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Abstract 

Much  of  the  recent  research  in  ocean  acoustics  has  focused  on  developing  meth¬ 
ods  to  exploit  the  effects  that  the  sea  surface  and  seaf  oor  have  on  acoustic  propaga¬ 
tion.  Many  of  those  methods  require  detailed  knowledge  of  the  acoustic  properties 
of  the  seaf  oor  and  the  sound  speed  profle  (SSP),  which  limits  their  applicability. 
The  range-frequency  waveguide  invariant  describes  striations  that  often  appear  in 
plots  of  acoustic  intensity  versus  range  and  frequency.  These  range-frequency  stri¬ 
ations  have  properties  that  depend  strongly  on  the  frequency  of  the  acoustic  source 
and  on  distance  between  the  acoustic  source  and  receiver,  but  that  depend  mildly  on 
the  SSP  and  seaf  oor  properties.  Because  of  this  dependence,  the  waveguide  invari¬ 
ant  can  be  utilized  for  applications  such  as  passive  and  active  sonar,  time-reversal 
mirrors,  and  array  processing,  even  when  the  SSP  or  the  seaf  oor  properties  are 
not  well  known.  This  thesis  develops  a  framework  for  understanding  and  calculat¬ 
ing  the  waveguide  invariant,  and  uses  that  framework  to  develop  signal  processing 
techniques  for  the  waveguide  invariant. 

A  method  for  passively  estimating  the  range  from  an  acoustic  source  to  a  re¬ 
ceiver  is  developed,  and  tested  on  experimental  data.  Heuristics  are  developed  to 
estimate  the  minimum  source  bandwidth  and  minimum  horizontal  aperture  required 
for  range  estimation. 

A  semi-analytic  formula  for  the  waveguide  invariant  is  derived  using  WKB  ap¬ 
proximation  along  with  a  normal  mode  description  of  the  acoustic  f  eld  in  a  range- 
independent  waveguide.  This  formula  is  applicable  to  waveguides  with  arbitrary 
SSPs,  and  reveals  precisely  how  the  SSP  and  the  seaf  oor  ref  ection  coeff  cient  af¬ 
fect  the  value  of  the  waveguide  invariant. 

Previous  research  has  shown  that  the  waveguide  invariant  range-frequency  stri¬ 
ations  can  be  observed  using  a  single  hydrophone  or  a  horizontal  line  array  (HLA) 
of  hydrophones.  This  thesis  shows  that  traditional  array  processing  techniques  are 
sometimes  inadequate  for  the  purpose  of  observing  range-frequency  striations  us¬ 
ing  a  HLA.  Array  processing  techniques  designed  specif  cally  for  observing  range- 
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frequency  striations  are  developed  and  demonstrated. 

Finally,  a  relationship  between  the  waveguide  invariant  and  wavenumber  inte¬ 
grations  is  derived,  which  may  be  useful  for  studying  range-frequency  striations  in 
elastic  environments  such  as  ice-covered  waveguides. 

Thesis  Supervisor:  Henrik  Schmidt 

Title:  Professor  of  Mechanical  and  Ocean  Engineering 
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Chapter  1 

Introduction 


Sound  provides  an  unparalleled  means  of  sensing  in  the  vast,  dark  oceans  that  cover 
70%  of  the  Earth.  Only  with  sound  can  one  detect  a  submerged  whale  from  hun¬ 
dreds  of  miles  away  [55],  or  measure  the  average  temperature  of  an  entire  ocean 
in  a  matter  of  hours  [65].  Although  sound  is  a  powerful  method  for  sensing  in 
the  ocean,  it’s  an  indirect  one:  the  desired  information  must  be  extracted  from  an 
acoustic  signal. 

Humans  do  this  type  of  acoustic  sensing  on  a  daily  basis  (albeit,  above  the  wa¬ 
ter).  A  person  can  determine  the  type  of  building  he  is  walking  in  —  a  large  stone 
cathedral,  or  an  off  ce  building?  —  based  solely  on  the  sound  of  his  footsteps.  He 
could  even  guess  roughly  how  large  the  cathedral  is,  based  on  how  much  rever¬ 
beration  he  hears.  But  to  answer  questions  like  these  quantitatively  —  how  many 
meters  wide  is  the  cathedral?  —  quantitative  methods  must  be  used. 

It  is  usually  not  obvious  how  to  extract  the  desired  information  from  the  raw 
acoustic  signal,  which  is  nothing  more  than  a  time-series  of  pressure  recorded  by 
a  hydrophone.  The  signal  usually  contains  not  only  the  information  of  interest, 
but  also  other  information  that  is  not  of  interest.  The  sound  caused  by  crashing 
surface  waves  may  not  be  of  interest  when  one  is  trying  to  detect  the  presence 
of  whale  vocalizations.  Conversely,  the  presence  of  whale  vocalizations  may  not 
be  of  interest  when  one  is  trying  to  detect  crashing  surface  waves.  One  of  the 
main  challenges  in  ocean-acoustic  sensing  is  to  create  signal  processing  techniques 
that  distill  out  the  information  of  interest  (e.g.,  whale  vocalizations)  while  ignoring 
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everything  else  (e.g.,  sound  of  waves  crashing). 

How  does  one  go  about  creating  such  signal  processing  techniques?  Some¬ 
times  this  can  be  done  with  standard  engineering  techniques,  such  as  using  a  low- 
pass  f  Iter  to  reject  high-frequency  noise.  But  often  in  ocean  acoustics,  the  signal 
processing  techniques  come  from  an  understanding  of  the  acoustics  itself. 

For  example,  ray  theory  reveals  that  some  of  the  acoustic  energy  in  the  ocean 
propagates  without  ever  being  affected  by  the  seaf  oor.  If  one  was  trying  to  acous¬ 
tically  determine  the  sound  speed  prof  le  (SSP)  of  the  water  column,  and  the  prop¬ 
erties  of  the  seaf  oor  were  unknown  and  not  of  interest,  one  could  use  ray  theory 
develop  a  signal  processing  technique  that  extracts  only  the  features  of  the  sig¬ 
nal  that  were  not  affected  by  the  seaf  oor.  The  output  of  such  a  signal  processing 
technique  could  then  be  used  to  infer  properties  of  the  water  column  without  ever 
having  to  be  concerned  with  the  properties  of  the  seaf  oor. 

This  thesis  focuses  on  a  feature  of  acoustic  propagation  in  a  waveguide  known 
as  the  “waveguide  invariant.”  The  waveguide  invariant  is  a  parameter  denoted  by 
/3  which  summarizes  the  frequency  dependence  of  a  waveguide’s  Green’s  function. 
The  waveguide  invariant  is  important  not  only  because  it’s  a  fundamental  property 
of  the  ocean  acoustic  waveguide,  but  also  because  it  can  be  exploited  to  extract  use¬ 
ful  information  from  an  acoustic  signal  when  the  details  of  the  ocean  environment 
(the  SSP  and  seaf  oor  properties)  are  unknown. 

Many  signal  processing  technique  that  utilize  the  physics  of  waveguide  propa¬ 
gation,  such  as  matched  f  eld  processing,  require  detailed  knowledge  of  the  ocean 
environment.  The  accuracy  of  such  methods  tends  to  be  very  sensitive  to  mis¬ 
matches  between  the  actual  environment  and  the  assumed  environment.  In  contrast, 
signal  processing  techniques  based  on  the  waveguide  invariant  often  require  only  a 
minimal  amount  of  knowledge  about  the  environment.  Because  of  this,  techniques 
based  on  the  waveguide  invariant  are  applicable  even  when  the  details  of  the  SSP 
or  the  seaf  oor  are  unknown. 

The  waveguide  invariant  has  been  used  for  a  wide  range  of  applications  such 
as:  passive  range  estimation  [58,  56,  62],  matched  f  eld  processing  [59,  23],  active 
sonar  [49,  27,  25],  array  processing  [37,  57,  66],  time-reversal  mirrors  [32,  54,  38, 
39],  and  more  (See  Appendix  A  for  a  complete  review  of  the  waveguide  invari- 
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ant  literature).  Despite  its  many  uses,  the  waveguide  invariant  is  not  well  studied 
compared  to  other  ocean-acoustic  phenomena. 

The  goal  of  this  thesis  is  to  improve  upon  the  latest  understanding  of  the  wave¬ 
guide  invariant  and  related  signal  processing  so  that  future  applications  of  the  wave¬ 
guide  invariant  can  be  better  executed. 

1.1  Thesis  outline  and  contributions 

Chapters  3  -  6  are  the  original  research  contributions  of  this  thesis. 

Chapter  2  This  chapter  provides  the  technical  background  for  the  thesis;  it  ex¬ 
plains  normal  mode  propagation  in  a  range-independent  waveguide,  and  pro¬ 
vides  a  detailed  introduction  to  the  range-frequency  waveguide  invariant. 

Chapter  3  This  chapter  presents  one  of  the  most  straight-forward  applications 
of  the  waveguide  invariant:  estimating  the  range  to  a  broadband  acoustic 
source  in  a  shallow-water  waveguide  using  a  single  acoustic  receiver  towed 
directly  toward  the  acoustic  source.  Previous  research  has  shown  that  a  two- 
dimensional  Fourier  transform  can  be  used  to  extract  information  about  the 
slope  of  the  waveguide  invariant  striations.  This  chapter  extends  that  research 
by  relating  the  ocean-acoustic  environmental  parameters  to  the  signal  pro¬ 
cessing  parameters,  and  applying  the  resulting  signal  processing  technique 
to  estimate  the  range  to  the  acoustic  source  (assuming  /3  =  1).  Heuris¬ 
tics  are  developed  to  estimate  the  minimum  source  bandwidth  and  minimum 
horizontal  aperture  required  for  range  estimation.  The  range  estimation  al¬ 
gorithm  is  tested  on  experimental  and  simulated  data  for  source  ranges  of 
500-2200  m  and  frequencies  from  350  to  700  Hz.  The  algorithm  is  accurate 
to  within  approximately  25%  for  the  cases  tested  and  requires  only  a  minimal 
amount  of  a  priori  knowledge  about  the  environment. 

Chapter  4  This  chapter  addresses  the  question:  “Why  is  the  waveguide  invari¬ 
ant,  invariant?”  Although  it  has  been  shown  analytically  that  (3  «  1  for  ideal 
waveguides,  numerical  and  experimental  results  have  revealed  that  /3  ~  1 
for  many  (but  not  all)  realistic  shallow-water  waveguides  as  well.  Ocean- 
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acoustic  techniques  that  utilize  the  waveguide  invariant,  such  as  those  pre¬ 
sented  in  Ch.  3,  often  assume  (3  =  1.  Therefore,  it  is  important  to  under¬ 
stand  why  (3  «  1  in  many  realistic  waveguides,  and  when  f3  deviates  greatly 
from  1. 

This  chapter  presents  a  method  for  calculating  (3  using  a  modal  WKB  de¬ 
scription  of  the  acoustic  f  eld  in  a  range-independent  waveguide,  which  re¬ 
veals  a  straightforward  relationship  between  the  SSP  and  (3.  That  relationship 
is  used  to  illustrate  why  non-uniformities  in  the  SSP  sometimes  have  such  a 
small  effect  on  f3  and  under  what  circumstances  the  non-uniformities  will 
have  a  large  effect  on  /3.  The  method  relies  on  implicit  differentiation  and 
thus  does  not  explicitly  solve  for  the  horizontal  wavenumbers  of  the  modes, 
making  it  applicable  to  waveguides  with  arbitrary  sound  speed  prof  les  and 
f  uid  bottom  half-spaces.  Several  examples  are  given,  including  an  analytic 
estimate  of  (3  in  a  Pekeris  waveguide. 

Chapter  5  This  chapter  generalizes  a  previously  known  result  which  showed  that 
a  horizontal  line  array  (HLA)  can  be  used  to  observe  waveguide  invariant  stri- 
ations  from  one  source  while  rejecting  noise  from  other  sources.  It  is  shown 
that  array  weights  commonly  used  for  planewave  beamforming,  such  as  uni¬ 
form  weights,  can  have  the  unintended  effect  of  suppressing  some  the  desired 
striations.  Insights  gained  from  Ch.  4  are  used  to  show  that  even  when  noise 
is  not  present,  array  processing  can  still  be  useful  because  it  can  suppress 
components  of  the  striation  pattern  that  are  not  useful  for  acoustic  sensing. 
Experimental  data  is  used  to  illustrate  the  ability  of  an  array  processor  to 
preserve  the  desired  striations  while  rejecting  noise. 

Chapter  6  This  chapter  shows  that  although  the  waveguide  invariant  is  typically 
def  ned  in  terms  of  normal  modes  or  ray  theory,  it  can  also  be  related  to  the 
wavenumber-integration  method  for  calculating  the  acoustic  f  eld  in  a  waveg¬ 
uide.  The  Wiener-Khinchin  Theorem  is  used  to  show  that  the  autocorrelation 
of  the  wavenumber-integration  kernel,  when  plotted  versus  wavenumber  dif¬ 
ference  and  frequency,  contains  striations  that  can  be  described  by  the  waveg¬ 
uide  invariant. 
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Chapter  7  This  chapter  summarizes  and  concludes  the  thesis.  Several  specif  c 
suggestions  for  future  areas  of  research  are  given. 

Appendix  A  This  appendix  is  an  extensive  but  brief  review  of  the  waveguide 
invariant  literature,  most  of  which  is  not  required  to  support  the  conclusions 
of  the  thesis  but  is  included  for  the  interested  reader. 


17 


18 


Chapter  2 

Background 


This  chapter  brief  y  reviews  the  aspects  of  ocean  acoustics  that  are  relevant  for  this 
thesis  —  normal  modes  in  a  range-independent  waveguide  —  and  then  provides  a 
detailed  introduction  to  the  range-frequency  waveguide  invariant. 


2.1  Acoustic  propagation  in  a  waveguide 

We  restrict  our  attention  the  acoustic  f  eld  caused  by  a  point  source  in  a  range- 
independent  ocean  waveguide.  As  far  the  acoustics  is  concerned,  the  environment 
is  fully  described  by  the  sound  speed  prof  le  (SSP)  in  the  water  column,  and  the 
sound  speed  and  density  prof  le  of  the  seaf  oor.  If  one  is  only  interested  in  the 
acoustic  f  eld  in  the  water  column,  then  the  seaf  oor  can  be  fully  characterized  by 
its  ref  ection  coeff  cient  as  function  of  horizontal  wavenumber  and  frequency. 

As  simple  as  that  physical  model  may  seem,  a  tremendous  amount  of  research 
was  required  to  develop  techniques  to  accurately  calculate  the  acoustic  f  eld  in  such 
a  waveguide.  The  acoustic  f  eld  be  calculated  and  conceptually  understood  using 
ray  theory,  normal  modes,  or  wavenumber  integration.  Other  methods  can  be  used 
for  calculating  the  acoustic  f  eld,  such  as  f  nite  element  analysis  or  the  parabolic 
equation,  but  those  methods  do  not  lend  themselves  to  understanding  the  underlying 
acoustics  of  propagation  in  a  waveguide. 

The  concepts  in  this  thesis  are  most  readily  understood  using  the  normal  mode 
description  of  the  acoustic  f  eld. 
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2.2  Normal  modes 


This  section  explains  the  aspects  of  normal  modes  that  are  important  for  under¬ 
standing  the  waveguide  invariant,  including  some  obscure  topics  not  typically  cov¬ 
ered  in  textbooks.  A  full  derivation  and  interpretation  of  normal  modes  can  be 
found  in  [33,  Chs.  2  and  5],  as  well  as  [60]  and  [10]. 

When  the  acoustic  source  and  receiver  are  separated  by  a  distance  greater  than 
a  few  water  depths,  the  complex  pressure  as  a  function  of  range  from  the  source  r, 
depth  z,  and  frequency  co  can  be  written  as  [33,  Eq.  5.14] 


p{r,  z,  u) 


M 


W4  H  MZsWmiz)  r— 

P\Zs)v  oTTV  m=l  V  Krm 


(2.1) 


where  zs  is  the  source  depth,  krrn  is  the  horizontal  wavenumber  of  mode  rn,  and 
ipmiz)  is  the  mode  function  of  mode  m.  krm  and  ^m{z)  usually  depend  on  oj,  the 
sound  speed  prof  le  c(z),  and  seaf  oor  properties,  and  are  obtained  by  solving  the 
depth-separated  wave  equation  [33,  Eq.  5.90] 


along  with  the  boundary  conditions  at  the  top  and  bottom  of  the  waveguide,  which 
depend  on  the  acoustic  properties  of  the  sea  surface  and  seaf  oor.  Equation  2.2  is  a 
second-order  differential  equation,  for  which  a  discrete  number  of  solutions  exist, 
each  one  denoted  by  the  mode  number  m. 

krm  and  'Wm(z)  can  be  calculated  analytically  in  an  ideal  waveguide  [33,  Sec.  2.4.4] 
and  semi-analytically  in  a  Pekeris  waveguide  [33,  Sec.  2.4.5].  Semi-analytic  solu¬ 
tions  exist  for  other  situations  as  well,  such  as  surface-trapped  modes  in  an  n2- 
linear  waveguide  ([33,  Sec.  2.5.1],  [11,  Sec.  6.6.1]),  although  the  solutions  are  in 
terms  of  non-elementary  math  functions.  (See  [9,  Sec.  3]  and  [10,  Sec.  4.4.1  ]  for  a 
complete  list  of  known  analytic  solutions.) 

For  waveguides  with  complicated  sound  speed  prof  les  and  non- vacuum  seaf  oors 
(i.e.,  fuid  or  solid),  krm  and  vpn(z)  do  not  have  exact  analytic  solutions.  They  can 
be  calculated  numerically  using  ocean  acoustic  software  such  as  Kraken  [48].  Or 
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they  can  be  calculated  approximately  using  the  WKB  approximation  or  perturba¬ 
tion  theory  ([50,  Sec.  l.A],  [36,  Sec.  3.1]).  The  WKB  approximation  provides 
insight  into  how  the  environment  (SSP,  seafoor)  affects  krm  and  V'm(z),  and  is 
discussed  in  Sec.  4.2  and  Appendix  B. 

2.2.1  Basic  interpretation  of  normal  modes 

Regardless  of  how  the  krm  values  are  calculated,  the  complex  pressure  f  eld  in 
Eq.  (2. 1)  is  a  sum  of  modes,  each  one  propagating  with  its  own  horizontal  wavenum¬ 
ber  krm.  In  this  thesis  we  assume  the  acoustic  medium  has  no  attenuation  so  that  all 
the  krm  are  real,  though  many  of  the  concepts  presented  could  likely  be  extended 
to  the  case  of  mildly  attenuating  media. 

Small  changes  in  the  SSP  or  the  bottom  boundary  conditions  (seaf  oor  proper¬ 
ties)  will  usually  lead  to  small  changes  in  the  values  of  the  horizontal  wavenumbers 
and  the  shapes  of  mode  functions.  But  because  krm  is  multiplied  by  the  range  r 
in  Eq.  (2.1),  which  is  typically  on  the  order  of  1000  m,  small  changes  in  the  krm 
lead  to  large  changes  in  the  structure  of  the  complex  pressure  f  eld.  An  example  of 
this  is  illustrated  in  Figs.  2.4  (b)  and  (e).  For  this  reason,  inversion  schemes  such 
as  matched  f  eld  processing  that  use  the  complex  pressure  level  as  measured  at  a 
specif  c  range,  depth,  and  frequency,  are  extremely  sensitive  to  environmental  mis¬ 
match.  Small  modeling  errors  (be  it  an  incorrect  environmental  parameter  value 
like  the  speed  of  sound  in  a  sediment  layer,  or  an  incorrect  parametrization  of  the 
environment  like  not  modeling  enough  sediment  layers)  can  prevent  an  inversion 
scheme  from  working  correctly. 

Another  quantity  of  interest  is  the  scalar  acoustic  intensity,  which  is  the  square 
of  the  magnitude  of  the  complex  pressure.  The  acoustic  intensity  as  a  function  of 
range  can  be  thought  of  as  an  interference  pattern  between  all  of  the  modes  [33, 
Sec.  2.4.4-5],  and  will  be  discussed  further  in  Sec.  2.3. 

2.2.2  Types  of  modes  (terminology) 

Understanding  the  relationship  between  a  mode’s  horizontal  wavenumber  and  the 
SSP  is  important  for  understanding  the  waveguide  invariant.  A  brief  summary  of 
the  relationship  is  given  here,  but  more  detailed  descriptions  can  be  found  in  [60, 
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Sec.  2.9],  [33,  Secs.  2.5.2  and  5.6],  and  [11,  Sec.  6.7],  The  terminology  used  here 
is  from  [33,  Sec.  1.4.1]  where  it  was  used  in  the  context  of  acoustic  rays. 

Surface-Reflected  Bottom-Reflected  ( SRBR )  modes  are  modes  that  have 
— -1 - —  <  krrn  <  -fr~  where  cmax  is  the  maximum  sound  speed  in  the  watercolumn 
(not  including  the  seaf  oor),  and  cseafoor  is  the  sound  speed  of  the  seafoor.  (For  a 
vacuum  or  rigid  seaf  oor,  use  cseaf  oor  =  oo). 

The  term  non-SRBR  will  be  used  to  refer  to  any  mode  that  is  not  an  SRBR 
mode.  Non-SRBR  modes  always  have  -fr-  <  krm  <  -f-  where  cmjn  is  the  mini- 

cmax  cmin 

mum  sound  speed  on  the  watercolumn.  Examples  of  non-SRBR  modes  are  surface- 
trapped  modes,  bottom-trapped  modes,  and  waterborne  modes. 


LJ  C 0  CO 


Csea  floor  Cmax  Cmin 

| < - SRBR - ►  - Non-SRBR ►  ] 

Figure  2.1:  Illustration  showing  the  relationship  between  the  sound  speed  prof  le, 
the  horizontal  wavenumber,  and  the  type  of  mode.  The  vertical  axis  is  depth  and 
the  horizontal  axis  is  wavenumber.  The  horizontal  axis  also  represents  the  hori¬ 
zontal  wavenumber  krm  of  any  given  mode.  Modes  with  horizontal  wavenumbers 
between  —J - —  and  -r-  are  surface-ref  ected  bottom-ref  ected  (SRBR)  modes. 

cseaf  oor  Cmax 


2.2.3  The  modal  sum 

As  shown  in  Eq.  (2.1),  the  total  acoustic  pressure  is  a  summation  of  terms  from 
M  normal  modes.  This  modal  sum  has  an  interesting  property  that  is  essential  for 
understanding  the  waveguide  invariant,  but  is  often  not  discussed. 

In  Chuprov’s  original  derivation  of  the  waveguide  invariant  [18],  he  states  that 
the  acoustic  f  eld  in  a  range-independent  ocean  is  typically  “composed  of  a  limited 
number  of  groups  of  modes  with  close  [mode]  numbers.”  This  subsection  explains 
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why  that  is  the  case,  by  summarizing  the  argument  presented  in  [20,  Sec.  3.1.6].1 

As  discussed  in  [33,  Sec.  5.2.1],  the  mode  number  m  is  a  monotonic  function  of 
kr  (this  property  can  also  be  understood  using  the  WKB  approximation  in  Ch.  4). 
Thus,  a  particular  range  of  kr  values  corresponds  to  a  group  of  adjacent  mode 
numbers.  The  quote  from  Chuprov  in  the  previous  paragraph  can  then  be  stated 
as:  the  complex  pressure  f  eld  at  a  given  range  is  usually  dominated  by  a  group 
of  modes  that  lie  within  a  small  range  of  kr  values  —  or  possibly  a  few  groups  of 
modes  with  a  few  distinct  ranges  of  kr  values. 

Each  term  of  the  summation  in  Eq.  (2.1)  represents  a  mode’s  contribution  to  the 
total  complex  pressure  at  a  f  xed  range,  depth,  and  frequency.  That  contribution  is 
a  complex  number  with  a  magnitude  and  phase,  which  can  be  plotted  as  a  vector  in 
the  complex  plane.  The  vector  corresponding  to  the  f  rst  mode  can  be  plotted  from 
the  origin  of  the  complex  plane.  The  vector  corresponding  to  the  second  mode  can 
plotted  from  the  end  of  the  f  rst  mode’s  vector;  the  third  one  plotted  from  the  end 
of  the  second  one,  and  so  on,  as  shown  in  Fig.  2.2  (a).  The  end  of  the  vector  from 
the  M-th  mode  is  then  the  total  complex  pressure  level  at  the  fxed  range,  depth, 
and  frequency. 

From  Eq.  (2.1),  it  can  be  seen  that  the  phase  difference  between  term  (m)  and 
term  (m  —  1)  is  r(krm  —  Av(m-i))-  Groups  of  modes  that  approximately  satisfy  the 
equation 

r{krm  -  Av(m_r))  =  2nN  (2.3) 

where  N  is  an  integer  will  all  have  approximately  the  same  phase  and  thus  will  add 
constructively  (they  will  all  be  pointing  in  the  same  direction  in  the  complex  plane); 
all  other  modes  will  be  distributed  in  “random”  directions  in  the  complex  plane 
and  will  approximately  cancel  each  other  out,  roughly  speaking.  This  argument  is 
similar  to  the  argument  used  to  justify  the  stationary  phase  method  of  integration. 

One  can  write  Eq.  (2.3)  as 


dkrm  _  2ttN 
dm  r 


(2.4) 


to  see  that  the  complex  pressure  level  is  dominated  by  a  few  groups  of  modes  (each 


*[18]  actually  referenced  [20,  Ch.  2],  but  that  was  likely  an  error  because  the  relevant  material 
appears  to  be  in  [20,  Ch.  3] 
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(a) 
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real  part 


(b) 


Figure  2.2:  (a)  Graphical  representation  of  the  mode  summation  for  a  f  xed  range, 
depth  and  frequency.  This  plot  shows  the  complex  pressure  in  the  deep-water 
waveguide  in  Fig.  2.6  at  a  range  of  80  km  and  100  Hz.  Each  line  segment  rep¬ 
resents  the  contribution  from  a  single  mode  number  (labeled),  beginning  with 
the  frst  modes  at  (0,0*)  and  ending  with  the  last  mode  about  (2.75,-3.25*). 
The  total  complex  acoustic  pressure  is  primarily  determined  by  modes  160  to 
170.  (b)  A  plot  of  the  phase  difference  between  terms  of  the  modal  summation: 
r(krm  —  fcr(m_1))/( 2ir).  The  phase  difference  for  modes  160  to  170  is  approx¬ 
imately  —2,  an  integer,  so  those  modes  all  have  approximately  the  same  phase 
(direction  in  the  complex  plane)  and  dominant  the  mode  sum. 
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one  corresponding  to  a  particular  value  of  N)  along  the  curve  of  kr  versus  m. 
Figure  2.2  (b)  is  a  plot  of  r(^^)/(27r)  versus  m  for  a  deep-water  waveguide  at 
80  km  in  range. 

The  concepts  in  the  present  subsection  are  illustrated  in  more  detail,  and  with 
more  mathematical  rigor  in  [20,  Sec.  3.1.6],  which  also  discusses  how  this  inter¬ 
pretation  of  the  modal  sum  relates  to  ray  theory. 

2.3  The  range-frequency  waveguide  invariant 

Historical  context  and  philosophical  note 

The  Russian  literature  on  ocean  acoustics  is  noticeably  different  than  the  Western 
(i.e.,  American  and  European)  literature.  Each  group  —  Russian  and  Western  — 
has  a  slightly  different  way  of  thinking  about  ocean  acoustics,  and  each  way  has  its 
own  strengths  and  weaknesses. 

One  of  the  lessons  learned  from  this  thesis  was  that  much  can  be  gained  by 
drawing  from  the  strengths  of  the  two  different  ways  of  thinking  about  ocean  acous¬ 
tics,  which  complement  each  other  very  well.  For  this  reason,  this  author  often  lists 
several  references  for  a  single  concept. 

The  waveguide  invariant,  the  topic  of  this  thesis,  was  discovered  by  a  Rus¬ 
sian  acoustician  named  S.D.  Chuprov  in  1982  [18].  Because  of  the  Cold  War 
between  the  United  States  and  Russia,  researchers  in  the  United  States  (and  pre¬ 
sumably  other  Western  countries)  did  not  know  about  the  waveguide  invariant  until 
the  1990s. 

Chuprov’s  original  paper  on  the  waveguide  invariant,  Ref.  [18],  has  been  trans¬ 
lated  to  English.2  Chuprov’s  paper  [18]  is  extremely  dense  —  covering  much  ma¬ 
terial  not  even  mentioned  in  this  thesis  —  and  requires  a  very  deep  understanding 
of  ocean  acoustics. 

2.3.1  Range-frequency  striations 

This  thesis  uses  the  term  “waveguide  invariant”  to  refer  to  the  range-frequency 
waveguide  invariant.  In  other  literature,  “waveguide  invariant”  sometimes  refers 

2This  author  obtained  a  copy  from  the  library  at  the  NATO  Undersea  Research  Centre. 


25 


to  a  slightly  more  general  concept  that  describes  not  only  changes  in  range  and 
frequency,  but  also  changes  environmental  parameters  such  as  waveguide  depth 
([18],  [34,  Sec.  2.4.6]). 

Sections  2.3.1  -  2.3.3  of  this  thesis  draw  heavily  from  [34,  Sec.  2.4.6]  and  [11, 
Sec.  6.7.2],  both  of  which  have  excellent  derivations  and  discussions  of  the  waveg¬ 
uide  invariant.  The  derivation  in  this  chapter  attempts  to  complement  those  two 
sources  by  presenting  the  concepts  in  a  slightly  different  manner,  and  by  explic¬ 
itly  showing  many  of  the  intermediate  steps  of  the  derivation.  This  section  uses 
normal  modes  to  describe  the  waveguide  invariant,  but  the  waveguide  invariant  can 
also  be  described  using  ray  theory  (see  Appendix  A.2).  Also,  Chapter  6  presents 
some  original  research  relating  the  waveguide  invariant  to  wavenumber  integration. 

We  start  out  with  an  empirical  observation:  A  2-d  plot  of  acoustic  intensity  in 
a  waveguide  as  a  function  of  range  and  frequency  (I(r,  cu)  =  \p(r,  cu)|2)  for  fxed 
source  and  receiver  depths  often  contains  striations.  This  plot  is  what  one  would 
obtain  if  he  or  she  plotted  a  spectrogram  from  an  acoustic  receiver  being  towed 
directly  away  from  a  broadband  acoustic  source  in  a  range  independent  waveguide. 

I(r,u )  versus  r  for  a  fxed  u  is  the  intermodal  interference  pattern  discussed 
in  [33,  Sec.  2.4.4-5],  so  the  striations  in  I(r,  u)  are  a  result  of  how  that  intermodal 
interference  pattern  changes  with  frequency.  The  striations  in  I(r,u )  will  be  re¬ 
ferred  to  as  waveguide  invariant  striations,  and  can  be  observed  under  a  wide  range 
of  environmental  conditions. 

A  simple  example  of  these  waveguide  invariant  striations  is  shown  in  Fig.  2.3  (a), 
which  shows  the  simulated  acoustic  intensity  in  a  Pekeris  waveguide  (source  at 
zs  =  40  m,  receiver  at  z  =  20  m,  cwater  =  1500  m/s,  bottom  =  1700  m/s, 
Pbottom  =  1750  kg/m3  —  but  all  Pekeris  waveguides  will  have  similar  striations). 
More  complicated  examples  are  shown  in  Figures  2.4  (c)  and  (d),  which  show  the 
waveguide  invariant  striations  for  two  shallow-water  waveguides  with  non-uniform 
sound  speed  prof  les.  Figure  2.5  shows  striations  from  experimental  data  collected 
in  a  shallow-water  waveguide  (see  Ch.  3  for  details).  Figure  2.6  shows  simulated 
waveguide  invariant  striations  for  a  deep-water  waveguide  with  a  Munk  sound 
speed  prof  le.  These  striations  from  a  deep-water  waveguide  are  different  look¬ 
ing  than  those  from  shallow-water  waveguides  shown  in  Figs.  2.3  and  2.4,  and  will 
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be  discussed  later. 


These  plots  show  that  similar  striation  patterns  appear  in  many  different  shallow- 
water  range-independent  environments  (with  the  deep-water  waveguide  having  a 
different  pattern).  In  shallow  water,  the  striation  pattern  is  often  not  strongly  af¬ 
fected  by  the  details  of  the  sound  speed  prof  le  or  the  seaf  oor  properties:  thus  the 
term  waveguide  invariant. 


The  acoustic  pressure  in  a  range-independent  waveguide  can  be  written  as  [see 
Eq.  (2.1)] 

^ ^  piformT 

(2.5) 


p(r,  uj)  oc 


keeping  in  mind  that  krrn  and  'ipm{z)  depend  on  uj.  Note  that  although  the  pressure 
depends  on  z,  as  well  as  r  and  uj,  the  z  dependence  is  not  explicitly  written  on 
the  left-hand-side  of  Eq.  (2.5)  because  the  z  dependence  is  not  important  for  this 
derivation.  Defne: 


A 


m 


'Ipm  (2) 


1 

y/krrnT 


(2.6) 


The  pressure  can  then  be  written  as 


p(r,u)  oc^Ameikrmr 

m 


(2.7) 


Am  is  a  slowly  varying  function  of  r  and  c 0  compared  to  how  quickly  the  exponen¬ 
tial  term  in  Eq.  (2.5)  varies  with  r  and  uj,  and  so  Am’s  dependence  on  r  and  uj  will 
be  ignored.  Am  can  be  thought  of  as  the  local  amplitude  in  the  (r,  to)  plane. 


The  scalar  acoustic  intensity  is  the  pressure  [Eq.  2.7]  multiplied  by  its  complex 
conjugate,  which  can  be  written  as  3 


3  Assuming  Am  and  Ai  are  real.  If  Am  and  Ai  are  complex,  then  the  AmAi  in  the  result  must  be 
replaced  with  |  AmA1  +  AiA*m  |,  which  is  approximately  equal  to  AmAi  when  the  real  parts  of  Am 
and  Ai  are  much  larger  than  the  imaginary  parts. 
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Figure  2.3:  (a)  Plot  illustrating  waveguide  invariant  striations  in  a  simulated  Pekeris 
waveguide.  The  quantity  plotted  is  acoustic  intensity  (dB  with  an  arbitrary  refer¬ 
ence)  versus  range  and  frequency.  The  striations’  slopes  (pointing  toward  the  lower 
left)  are  described  by  the  waveguide  invariant  3.  (b)  Same  as  (a),  but  overlaid  with 
lines  corresponding  to  3  =  1/2.  (c)  same  as  (b)  but  with  lines  corresponding  to 
3  =  1.  (d)  same  as  (b)  but  with  lines  corresponding  to  3  =  3/2.  Note  that  the 
striations  are  best  described  with  3  =  1. 
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Figure  2.4:  Plots  illustrating  the  waveguide  invariant  for  two  different  simulated 
shallow-water  waveguides,  with  the  properties  of  one  waveguide  plotted  in  the  left 
column  and  the  properties  of  the  other  waveguide  in  the  right  column.  The  top  row 
shows  the  sound  speed  prof  les  (including  the  seaf  oor).  The  second  row  shows  the 
acoustic  intensity  (dB  with  an  arbitrary  reference)  versus  range  at  350  Flz,  with  the 
source  at  37  m  deep  and  the  receiver  at  10  m  deep.  The  plots  in  the  second  row  illus¬ 
trate  that  two  waveguides  with  similar  (but  not  identical)  sound  speed  prof  les  can 
have  dramatically  different  acoustic  f  elds.  The  plots  in  the  third  row  are  acous¬ 
tic  intensity  (dB  with  an  arbitrary  reference)  versus  range  and  frequency.  Both 
waveguides,  despite  their  differing  sound  speed  prof  les,  contain  similar-looking 
striations,  the  slopes  of  which  are  described  by  the  waveguide  invariant  j3. 
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Figure  2.5:  Acoustic  intensity  (dB  with  an  arbitrary  reference)  versus  range  and 
frequency  illustrating  waveguide  invariant  striations  with  experimental  data  from  a 
shallow-water  waveguide. 
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Figure  2.6:  Similar  to  Fig.  2.4,  but  for  a  deep-water  waveguide  with  a  Munk  sound 
speed  prof  le  and  plotted  with  different  range  and  frequency  scales.  The  source 
is  at  300  m  deep  and  the  receiver  is  3100  m  deep,  so  low  order  modes  contribute 
signifcantly  to  the  total  acoustic  intensity.  Note  that  the  striation  pattern  looks 
qualitatively  different  than  those  in  Figs.  2.3  (a)  and  2.4  (c)  and  (e). 
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(2.8) 


I{r,u)  =  \p(r,u)\2 


y,  y  (AmAi  cos(A kmir)) 

m  l 

(2.9) 

y  A2  +  y  AmAi  cos(Akmir) 

(2.10) 

q  m,l ;  m^l 


where  m  and  l  are  the  mode  number  indices,  A kmi  =  krm  —  kri  depends  on  oj, 
and  the  exponents  have  been  written  as  cosines  using  Euler’s  formula.  Inspection 
of  Eq.  (2.10)  reveals  that  the  acoustic  intensity  is  a  sum  of  cosines,  each  cosine 
resulting  from  the  interference  between  two  modes  and  having  the  amplitude  of  the 
product  of  those  two  modes’  amplitudes. 

Defne: 

Iml ( r ,  w)  =  AmAi  cos(A kmir)  (2.11) 

so  that  Imi(r ,  oS)  is  the  contribution  to  the  total  intensity  from  a  single  cosine  term. 
A  plot  of  a  Imi(r,uj)  versus  cu  and  r,  will  contain  striations,  as  shown  in  Fig.  2.7 
for  m  =  4  and  l  =  10. 

The  slope  of  those  striations  can  be  determined  by  f  nding  the  direction  in  the 
(r,  ui)  plane  for  which  the  intensity  does  not  change.  This  is  done  by  taking  a  f  rst- 
order  two-dimensional  Taylor  series  expansion  of  Imi(r,  u)  about  some  point  in  the 
(r,  u)  plane,  and  setting  the  result  equal  to  zero: 


dlmi(r,u ) 

- a - dr  + 

or 


dlmi(r,u ) 

- - - bu  =  0 

oui 


(2.12) 


Suj 

Then  the  slope  —  (the  direction  in  which  /,,,/  (r,  to)  does  not  change  value)  can  be 
dr 

solved  for: 


8lo 

Sr 


dlmi(r,u )  / dlmi(r,u ) 


dr 


du 


(2.13) 


From  the  perspective  of  a  person  walking  around  on  the  bumpy  surface  of 
Iml(r,uj),  Eqs.  (2.12)  and  (2.13)  can  be  stated  in  words  as:  If  I  am  at  a  location 
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Figure  2.7 :  A  plot  of  a  single  cosine  term  of  the  acoustic  intensity  for  the  waveguide 
shown  in  Fig.  2.4  (d),  computed  using  modes  4  and  10. 


(r,  u)  and  I  move  in  the  r  direction  by  some  small  amount  5r,  then  how  far  must  I 
move  in  the  u  direction  (6uj)  so  that  Imi  has  the  same  value  as  where  I  started?  A 
“striation”  is  a  line  or  curve  in  the  (r,  co)  plane  for  which  /,,,;  does  not  change. 
Inserting  Eq.  (2.11)  into  Eq.  (2.13)  yields 

5lo  _  AmAi  •  (Akmi(uj))  sin(Akmi(uj)  r)  _  (A kmi(u)) 

5r  rAmAi  ■  {dAk^(uj))  sin(A kmi(u)  r )  f  dAkmi(uj)\ 

V  duj  ) 


We  then  defne  the  waveguide  invariant  parameter  for  these  two  modes,  fimu 


as: 


1  A/rm/(cu) 

UJ  dAkmi{u)/du 


(2.15) 


The  reason  for  this  def  nition  is  not  immediately  obvious,  but  it  turns  out  to  be  a 
useful  one.  Using  Eq.  (2.15)  and  the  right  hand  side  of  Eq.  (2.14),  one  can  write 


~F~  =  fiml  ~~  ■  (2.16) 

dr  r 

So  the  slope  ^  of  a  striation  of  a  single  Imi  term  is  3rni  times  the  frequency  u 
divided  by  the  range  r. 

But  the  acoustic  intensity  I(r,ui)  consists  of  a  sum  of  Imi  terms,  and  each 
(5mi  only  describes  striations  in  a  single  Imi  term.  If  all  of  the  Irni  terms  have 
8mi  values  that  are  approximately  the  same,  then  each  term  will  have  striations 


33 


with  approximately  the  same  slopes,  and  so  the  total  intensity  I(r,ui)  will  have 
striations  described  by  a  scalar  value  j3  that  is  approximately  independent  of  m  and 
l.  Section  2.3.2  shows  that  in  an  ideal  waveguide,  j3mi  is  approximately  independent 
of  m,  l,  and  uj,  for  modes  far  from  cut-off. 

In  more  complicated  waveguides,  however,  j3mi  is  not  independent  of  m,  l,  and 
uj.  In  that  case,  a  more  sophisticated  view  of  the  waveguide  invariant  is  required  in 
order  to  understand  why  there  are  striations  in  I(r,uj)  even  though  I(r,u)  consists 
of  a  sum  of  Imi  terms,  each  of  which  has  a  different  / 3mi  value  (and  thus  different 
striations  slopes).  This  more  sophisticated  view  is  described  in  Sec.  2.3.3.  But  f  rst, 
Sec.  2.3.2  shows  that  l3mi  «  1  in  an  ideal  waveguide. 

Note  that  in  Eq.  2.16,  one  can  replace  uj  (the  temporal  frequency  in  radians  per 
second)  with  /  (the  temporal  frequency  in  cycles  per  second)  because  a  factor  of 
2-7t  will  appear  on  both  sides  of  the  equation  and  will  thus  cancel  out.  Also  note 
that  Eq.  2.16  is  approximate  because  the  dependence  of  Arn  on  r  and  uj  has  been 
ignored. 


If  i  <  1  (which  means  the  vertical  wavenumbers  are  much  smaller  than  the  total 
wavenumber,  or  equivalently  that  the  angle  of  propagation  is  close  to  horizontal), 
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then  we  can  Taylor  series  expand  y/1  —  x  to 


VT^x  =  1  -  |  +  0(x2).  (2.20) 

Inserting  x  in  Eq.  (2.20)  and  ignoring  ()(x2)  terms,  we  obtain 

krm~  fc(1_  2  (fed)  )  (221) 

where  we  used  the  approximate  sign  because  we  are  ignoring  all  terms  of  0(x2)  or 
higher.  Under  that  assumption  (and  replacing  k  with  -), 

Km  -  Kl «  ^  Q)  (m2  -  )2)  (2.22) 

This  can  be  inserted  into  Eq.  (2.15)  to  reveal  that  f3  «  1  in  this  case,  regardless 
of  the  values  of  rn,  l,  and  u  (bearing  in  mind  that  this  is  only  valid  for  modes 
with  small  vertical  wavenumbers  compared  to  their  horizontal  wavenumbers).  An 
alternative  derivation  expresses  the  same  result  in  terms  of  the  modal  propagation 
angles  [18,  Eq.  9]  [34,  Sec.  2.4.6]. 

Modes  far  from  cut-off  in  a  Pekeris  waveguide  have  horizontal  wavenumber 
differences  and  frequency  dependencies  that  are  similar  to  modes  in  an  ideal  waveg¬ 
uide.  To  see  this,  consider  a  plot  of  the  characteristic  equation  for  a  Pekeris  waveg¬ 
uide  [33,  Eq.  2. 167]  versus  horizontal  wavenumber,  and  how  that  plot  would  change 
with  frequency.  Because  of  this,  modes  far  from  cut-off  in  a  Pekeris  waveguide  also 
have  Pmi  «  1. 

Experimental  and  numerical  observations  suggest  that  3mi  ~  1  in  many  shallow- 
water  waveguides,  regardless  of  their  sound  speed  prof  les.  Qualitative  reasons  for 
this  are  suggested  in  [18];  Chapter  4  investigates  this  in  more  depth. 

2.3.3  The  waveguide  invariant  non-ideal  waveguides 

One  way  to  understand  the  waveguide  invariant  in  non-ideal  waveguides  is  to  ma¬ 
nipulate  Eq.  (2. 15)  to  def  ne  f3mi  in  terms  of  the  group  slowness  and  phase  slowness 
(the  reciprocal  of  group  speed  and  phase  speed,  respectively).  Group  and  phase 
speeds  are  discussed  in  [33,  Sec.  2.4.4]  and  [10,  Sec.  4.5]  (with  rather  different 
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points-of-view). 

The  phase  speed  for  a  mode  m  is 


u 

vm  ~  , 

A,  7*777, 

(2.23) 

and  the  group  speed  is 

du 

(2.24) 

u m  =  a, 

ukrm 

The  phase  slowness  is 

S  _  J_ 

^ p,m  — 

(2.25) 

Vm 

and  the  group  slowness  is 

s  _  J_ 

^ g,m  — 

(2.26) 

Comparing  those  defnitions  with  the  right-hand-side  of  Eq.  (2.15),  one  can  write: 


fiml  — 


AS p,ml  (^0 


(2.27) 


where  A Sp>mi  =  SPtm  —  Spj,  and  analogously  for  ASgmi.  (The  last  several  equa¬ 
tions  are  contained  in  original  waveguide  invariant  paper.  Ref.  [18],  but  the  notation 
and  interpretation  used  in  this  thesis  come  from  [34,  Sec.  2.4.6]  and  [66].) 

The  meaning  of  Eq.  (2.27)  can  be  illustrated  by  a  plot  of  Sg  versus  Sp  for  all 
the  1  through  M  modes  (at  a  f  xed  w),  as  is  done  by  Ref.  [18]  and  [34,  Sec.  2.4.6], 
This  plot  of  Sg  versus  Sp  can  be  thought  of  as  a  functional  relationship:  Sg  as  a 
function  of  Sp  with  both  Sg  and  Sp  being  parameterized  by  the  mode  number  m. 
This  function  will  be  denoted  as  Sg(Sp).  In  fact,  under  the  WKB  approximation, 
it  can  be  shown  that  Sg  is  a  function  of  Sp,  independent  of  u,  for  modes  that  do 
not  interact  with  the  seaf  oor.  This  is  mentioned  in  Ref.  [18]  and  can  be  seen  in 
Eq.  (4.36)  of  this  thesis. 

If  we  consider  only  regions  of  Sg{Sp)  where  Sg(Sp)  is  approximately  a  straight 
line  (mathematically  speaking,  take  a  f  rst  order  Taylor  series  expansion  of  Sg(Sp) 
about  some  f  xed  Sp  value),  then  3  is  def  ned  in  terms  of  the  slope  of  that  line: 


1  _  dSg 

J(Sp)=~dSP' 


(2.28) 
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Figure  2.8:  The  plots  in  the  left  column  correspond  to  waveguide  in  the  left  column 
in  Fig.  2.4,  and  the  plots  in  the  right  column  correspond  to  waveguide  in  the  right 
column  in  Fig.  2.4.  The  top  row  is  a  plot  of  group  speed  versus  phase  speed  at 
350  Hz,  with  each  dot  corresponding  to  a  mode.  The  bottom  row  is  a  plot  of  group 
slowness  versus  phase  slowness  at  350  Hz.  The  slope  of  a  line  connecting  points  of 
group  slowness  versus  phase  slowness  is  the  reciprocal  of  the  waveguide  invariant. 
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Figure  2.9:  Plots  of  (a)  group  speed  versus  phase  speed  and  (b)  group  slowness  ver¬ 
sus  phase  slowness  for  the  deep-water  waveguide  shown  in  Fig.  2.6,  at  50  Flz.  Each 
dot  corresponds  to  one  mode.  The  slope  of  a  line  connecting  points  in  the  group 
slowness  versus  phase  slowness  plot  is  the  reciprocal  of  the  waveguide  invariant. 
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This  3  no  longer  has  an  ml  subscript  because  it  depends  only  on  the  phase  slowness 
Sp  in  Sg(Sp),  as  opposed  to  depending  on  the  phase  slownesses  of  two  distinct 
modes  as  does  Eq.  (2.27).  When  calculated  this  way,  (3  corresponds  to  two  modes 
with  adjacent  modes  numbers.  3  may  depend  on  to  as  well,  as  will  be  discussed  in 
Ch.  4. 

It  may  seem  strange  (or  even  useless)  to  defne  3  using  only  adjacent  modes, 
given  that  the  acoustic  intensity  contains  components  from  all  mode  pairs,  includ¬ 
ing  those  with  m  and  l  values  that  are  not  close  to  each  other  [see  Eq.  (2.9)]; 
each  Imi(r,u )  will  have  striations  according  to  its  /3m[ .  However,  as  discussed  in 
Sec.  2.2.3,  the  acoustic  f  eld  at  f  xed  range  and  frequency,  is  dominated  by  groups 
of  modes  with  adjacent  mode  numbers.  That  provides  justif  cation  for  def  ning  3 
using  only  modes  with  adjacent  mode  numbers  [18],  This  author  is  not  aware  of  any 
research  attempting  to  quantify  the  validity  of  such  a  justif  cation.  But  evidently, 
Eq.  (2.28)  often  provides  an  accurate  description  of  range-frequency  striations. 

Figure  2.8  (a)  and  (c)  shows  plots  of  phase  speed  and  group  speed  for  the  two 
shallow-water  waveguides  discussed  earlier  in  this  chapter.  All  modes  with  phase 
speeds  less  than  1520  m/s  —  the  speed  of  sound  at  the  top  of  the  water  column 
—  are  non-SRBR  mode.  An  analogous  statement  can  be  made  for  the  phase  slow¬ 
nesses:  all  modes  with  Sp  =  >  152q  m/s  «  6.57  s/m  are  non-SRBR  modes.  It 

can  be  seen  in  Fig.  2.8  that  the  slope  of  the  line  connecting  non-SRBR  modes  is  sig- 
nif  cantly  different  than  the  slope  of  the  line  connecting  the  SRBR  modes.  Thus,  the 
non-SRBR  modes  and  the  SRBR  modes  each  correspond  to  signif  cantly  different 
values  of  3-  However,  because  there  are  many  more  SRBR  modes  than  non-SRBR 
modes  in  these  cases,  the  striation  pattern  observed  in  I ( r ,  u)  corresponds  mostly 
to  the  value  of  3  for  the  SRBR  modes(/3  «  1). 

Figure  2.9  is  the  same  as  Fig.  2.8,  but  for  the  deep-water  waveguide  shown  in 
Fig.  2.6.  In  this  case,  a  signif  cant  fraction  of  the  modes  are  non-SRBR  modes 
(which  in  this  case  are  modes  with  phase  speeds  less  than  about  1550  m/s).  That 
partially  explains  why  the  striations  in  Fig.  2.6  (c)  do  not  look  the  same  as  those 
in  Fig.  2.9(c)  and  (f):  non-SRBR  modes  tend  to  have  different  values  of  3  than  the 
SRBR  modes  do. 

Furthermore,  the  slope  of  the  striations  in  Fig.  2.6  (c)  appear  to  vary  depending 
on  the  range  and  frequency,  which  makes  sense  given  the  argument  in  Sec.  2.2.3: 
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different  locations  in  Fig.  2.9  are  dominated  by  different  groups  of  modes,  which 
have  different  values  of  3.  This  author  is  not  aware  of  any  research  directly  inves¬ 
tigating  how  the  value  of  /3  observed  in  the  local  striations  in  I(r,u)  may  change 
with  range  and  frequency  (See  Sec.  7.3  for  a  suggestion  on  how  to  investigate  this). 


In  many  cases,  however,  the  range-frequency  striation  pattern  is  well  described  by  a 

single,  dominant  value  of  3  that  is  approximately  independent  of  range,  frequency, 

5u 

and  mode  numbers.  In  such  cases,  the  striation  slopes  —  can  often  be  described 

or 

with  the  equation 

(2.29) 

or  r 

which  does  not  include  33  dependence  on  the  mode  numbers  m  and  l,  and  the 
frequency  to  (or  alternatively  Sp  and  to).  Eq.  (2.29)  is  a  useful  approximation  in 
many  situations.  However,  as  pointed  out  by  Rouseff  and  Spindel  in  [51],  it’s 
important  to  keep  in  mind  that  Eq.  (2.29)  is  an  approximation,  and  that  /3  is  not 
always  equal  to  a  single  value. 

Eq.  (2.29)  defnes  the  slope  of  a  striation  (assuming  the  striation  is  perfectly 
well  described  by  j3),  and  can  be  integrated  to  yield  the  equation  for  the  striation 
path  itself  [21,  Eq.  20]: 

to  =  coo  (2.30) 


where  coo  and  ro  represent  the  constant  of  integration,  and  def  ne  a  point  in  the 
(r,  to)  plane  where  the  striation  will  go  through.  When  the  source  and/or  receiver 
are  moving,  the  r  and  to  in  Eq.  (2.30)  can  be  parameterized  by  time,  as  shown  in 
[21,  Eq.  20]  by  D’  Spain  et  al. 


2.4  Chapter  summary 

The  total  acoustic  intensity  is  a  summation  of  cosine  terms,  each  representing  the 
contribution  from  a  pair  of  modes  denoted  by  m  and  l.  Each  cosine  term  —  each 
mode  pair  —  has  its  own  / 3mi  value  that  may  depend  on  to.  The  total  intensity  is 
the  addition  of  all  those  cosine  terms. 

Under  some  circumstances,  the  total  acoustic  pressure  and  thus  the  total  acous- 
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tic  intensity  is  dominated  by  mode  pairs  that  all  have  similar  values  of  j3mi.  This 
leads  to  the  approximation  that  (5  is  a  scalar  parameter. 

Understanding  that  the  acoustic  intensity  is  a  summation  of  cosines,  each  of 
which  has  different  (but  possibly  similar)  (3mi  value,  is  useful  for  understanding 
plots  of  I(r,  u).  This  understanding  may  also  be  useful  for  applying  the  waveguide 
invariant  concept  to  other  situations  such  as  array  processing,  time-reversal  mirrors, 
etc. 


This  chapter  reviewed  the  background  material  necessary  for  understanding  the 
original  work  presented  in  the  rest  of  this  thesis,  all  of  which  uses  /3  described  by 
normal  modes  in  a  range-independent  waveguide.  However,  much  more  research 
has  been  done  on  the  waveguide  invariant,  including  relationships  to  ray  theory  and 
extensions  to  range-  and  azimuthally-dependent  waveguides.  Appendix  A  is  a  brief 
but  thorough  review  of  the  waveguide  invariant  literature. 
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Chapter  3 

Robust  passive  range  estimation 
using  the  waveguide  invariant 


The  previous  chapter  explained  that  the  waveguide  invariant,  f3,  describes  range- 
frequency  striations.  In  this  chapter,  original  work  is  presented  that  utilizes  range- 
frequency  striations  to  estimate  the  range  to  an  acoustic  source  by  assuming  that 
j3  =  l.1  Because  this  chapter  presents  experimental  data  that  has  /3  «  1,  this 
chapter  also  serves  as  a  motivation  for  Ch.  4  which  investigates  why  /3  «  1  in  most 
shallow- water  waveguides. 

As  with  the  rest  of  this  thesis,  mode  numbers  are  denoted  by  m  and  l.  In 
contrast  to  the  rest  of  the  thesis,  this  chapter  writes  the  temporal  frequency  as  / 
(cycles/s)  instead  of  to  (radians/s):  u  =  2irf. 


3.1  Introduction 


Most  research  on  acoustic  source  localization  in  the  ocean  has  focused  on  using  co¬ 
herent  signal  processing  techniques  such  as  Matched  Field  Processing  [3]  (MFP). 
While  MFP  works  well  in  theory  and  in  numerical  simulation,  it  is  often  not  ap¬ 
plicable  to  real-world  situations  because  it  requires  very  accurate  knowledge  of  the 

'This  chapter  is  based  on  “Robust  passive  range  estimation  using  the  waveguide  invariant”  by 
Kevin  L.  Cockrell  and  Henrik  Schmidt  [The  Journal  of  the  Acoustical  Society  of  America.  May, 
2010], 
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environment  (e.g.,  sound  speed  prof  le  and  acoustic  properties  of  the  sea  f  oor)  in 
order  to  correctly  localize  the  source. 

Other  source  localization  methods  have  been  developed  that  require  much  less 
a  priori  knowledge  about  the  environment  by  using  the  waveguide  invariant.  The 
waveguide  invariant  has  been  applied  to  estimating  the  range  to  acoustic  sources  in 
various  circumstances  using  a  variety  of  signal  processing  schemes  such  as: 

•  Estimating  the  range  to  a  f  xed  acoustic  source  from  measurements  taken  by 
a  vertical  hydrophone  array  by  analyzing  the  MFP  sidelobe  behavior  [59], 

•  Estimating  the  closest-point-of-approach  of  a  moving  source  to  a  single  f  xed 
hydrophone  using  a  Plough  transform-like  technique  [56]. 

•  Estimating  source  trajectories  based  on  range-frequency  striations  in  an  ar¬ 
ray’s  beamformed  output  [62], 

In  this  chapter,  a  technique  is  developed  to  estimate  the  range  to  a  f  xed  acoustic 
source  from  the  acoustic  intensity  as  measured  over  a  window  of  ranges  and  fre¬ 
quencies,  /(r,  /).  The  technique  is  tested  on  experimental  data  that  was  obtained 
from  an  acoustic  receiver  towed  by  an  Autonomous  Underwater  Vehicle  (AUV) 
heading  directly  toward  the  acoustic  source.  Previous  research  papers  have  used 
the  two-dimensional  discrete  Fourier  transform  (2D-DFT)  of  I(r,  /)  to  estimate 
the  value  of  the  waveguide  invariant  when  the  source  range  was  known  [51,  66], 
but  did  not  address  how  to  choose  the  signal  processing  parameters.  This  chapter 
extends  that  research  to  using  the  2D-DFT  of  I (r,  /)  to  perform  range  estimation, 
and  develops  guidelines  for  choosing  the  signal  processing  parameters’  values. 

The  objective  of  the  present  work  is  to  investigate  the  issues  related  to  the  signal 
processing  that  is  required  for  range  estimation  using  the  2D-DFT  of  I (r,  /)  in  the 
context  of  performing  the  estimation  autonomously  (i.e.,  without  requiring  human 
interpretation  of  any  images),  and  to  perform  the  range  estimation  on  simulated  and 
experimental  data.  The  main  results  are: 

•  The  minimum  bandwidth  of  the  acoustic  source  and  minimum  range  win¬ 
dow  of  acoustic  intensity  measurements  required  for  accurate  range  estima¬ 
tion  can  be  determined  from  a  modest  knowledge  of  the  acoustic  waveguide 
parameters  before  any  acoustic  measurements  are  made. 
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•  A  relationship  between  the  signal  processing  parameters  and  the  ocean-acoustic 
waveguide  parameters  can  be  used  to  reject  much  of  the  noise  present  in  ex¬ 
perimental  data. 

•  Range  estimation  can  be  performed  robustly,  requiring  very  little  a  priori 
environmental  knowledge  (at  least,  for  the  data  sets  analyzed  this  chapter). 

•  A  range  estimate  accuracy  of  approximately  25%  is  achieved  with  the  exper¬ 
imental  data  set  used  in  this  chapter. 

Section  3.2  discusses  the  waveguide  invariant  in  the  context  of  passive  range 
estimation.  In  Sec.  3.3,  the  range  estimation  algorithm  based  on  a  2D-DFT  is  dis¬ 
cussed.  Section  3.4  applies  the  algorithm  to  simulated  and  experimental  data.  And 
f  nally,  Sec.  3.5  gives  a  summary  and  conclusion  of  this  chapter. 

3.2  Brief  review  of  range-frequency  striations 

Section  2.3  explained  that  a  plot  of  acoustic  intensity  versus  range  and  frequency 
due  to  a  broadband  source  in  a  waveguide,  I(r,  /),  will  exhibit  striations  that  are 
described  by  the  waveguide  invariant  (3-  As  discussed  in  Sec.  2.3  the  value  of  B  is 
unique  for  each  mode  pair,  but  /3  can  sometimes  be  approximated  as  being  a  single 
value  independent  of  the  mode  numbers  rn  and  l,  and  the  frequency  ui.  Under  that 
approximation,  specif  c  value  of  (3  that  will  be  observed  in  I(r,  co)  depends  on  the 
sound  speed  profle,  the  seafoor  properties,  and  the  source  and  receiver  locations 
[51,21,49,  59]. 

Empirical  and  numerical  observations  suggest  that  for  mode  pairs  where  both 
modes  are  surface-ref  ecting  bottom-ref  ection  (SRBR)  modes,  (3  «  1  (this  is  inves¬ 
tigated  further  Ch.  4).  So  if  the  acoustic  intensity  is  dominated  by  SRBR  modes, 
the  (3  observed  in  the  striation  pattern  will  likely  be  close  to  1 .  Throughout  this 
chapter,  for  both  the  simulated  and  experimental  data,  we  assume  (3  =  1.  The 
consequences  of  assuming  (3  =  1  are  discussed  shortly. 

Eq.  (2.29)  can  be  used  for  range  estimation  by  rewriting  it  as 

Sr 

r  =  (3.1) 

Sj 
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which  allows  for  one  to  estimate  the  range  of  the  acoustic  source  if  one  measures 
the  slopes  of  the  striations  and  assumes  a  value  of  f3.  The  effect  of  assuming  an 
incorrect  value  of  (3  can  been  seen  in  Eq.  (3.1).  If  the  true  value  of  (3  is  tflruc  and 
the  assumed  value  is  /3assumed,  then  the  range  estimates  will  be  incorrect  by  a  factor 

Of  /^assumed  / /^true  • 

In  order  to  estimate  the  range  to  the  source,  one  begins  by  calculating  /(r,  /) 
for  some  range  of  values  of  r  and  /.  In  a  simulation  one  can  calculate  I  (r,  /)  in  the 
frequency  domain  using  acoustic  simulation  software.  In  an  experiment  one  must 
estimate  the  power  spectrum  of  a  hydrophone’s  time  series  at  several  ranges.  In 
practice,  I(r,  /)  will  likely  be  the  spectrogram  of  a  time  series  of  acoustic  pressure 
obtained  by  moving  acoustic  an  receiver  radially  towards  or  away  from  the  acoustic 
source,  as  is  done  in  Sec.  3.4.2. 

To  estimate  the  source’s  location,  one  then  must  determine  the  slopes  of  the 
striations  (or  curved  paths  of  the  striations,  if  (3  /  1)  in  I(r,  /).  Because  of  the 
visually  striking  relationship  between  the  striation  slope  and  the  source’s  location, 
a  person  looking  at  I(r, /)  can  estimate  the  source’s  location  rather  easily  (see 
Fig.  2.3.  However,  the  present  work  focuses  on  techniques  that  perform  the  range 
estimation  autonomously  (i.e.,  without  the  beneft  of  having  a  person  to  visually 
interpret  I(r,  /)  or  its  2D-DFT). 


3.3  Using  the  two-dimensional  discrete  Fourier  transform 
for  range  estimation 

3.3.1  Outline  of  the  2D-DFT  technique 

The  slope  of  the  striation,  at  a  particular  range-frequency  combination  (r,  /) 
can  be  inserted  into  Eq.  (3.1)  to  estimate  the  range  to  the  source.  The  technique 
described  in  this  section  to  determine  the  slope  of  a  striation  in  7(r,  /)  looks  at  a 
small  local  region  (a  “window”)  of  I(r,  / )  and  assumes  that  all  of  the  striations 
within  that  window  have  the  same  slope.  This  is  similar  to  what  is  done  in  [4,  51, 
66],  but  in  those  papers  the  ranges  were  much  larger  than  the  ranges  used  in  the 
present  analysis,  so  the  slopes  of  the  striations  did  not  change  quickly  with  range. 
Consequently,  those  papers  did  not  focus  much  attention  on  how  to  choose  the  size 
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of  the  window.  Because  of  the  short  ranges  used  in  the  present  analysis,  the  slopes 
change  quickly  with  range  and  so  much  care  must  be  given  to  the  choosing  the 
window  size  in  order  to  ensure  that  the  slopes  do  not  change  too  much  inside  of  the 
window.  This  issue  is  discussed  in  Sec.  3.3.2. 

In  addition  to  the  striations  from  the  source,  I(r,  /)  will  also  contain  noise 
which  can  be  partially  eliminated  by  f  ltering.  The  spatial  cutoff  frequencies  of  the 
f  Iter  are  discussed  in  Sec.  3.3.2. 

The  range  is  then  estimated  based  on  the  slope  of  the  striations  in  the  window. 
This  process  is  repeated  for  several  windows  located  on  a  grid  in  the  (r,  /)  plane. 
Each  window  on  the  grid  will  produce  one  range  estimate,  all  of  which  can  then  be 
averaged  obtain  a  single,  robust  estimate.  Note  that  the  vertical  axis  of  the  window 
is  frequency,  so  it  requires  that  the  source  be  broadband.  The  horizontal  axis  of  the 
window  is  range,  and  so  it  requires  that  the  acoustic  f  eld  is  measured  along  a  line 
emanating  radially  from  the  acoustic  source. 

3.3.2  Determining  local  striation  angle  with  the  2D-DFT 

Denote  a  rectangular  window  of  I(r,  /)  bounded  by  (rmjn  <  r  <  rmax)  and 
(/min  <  /  <  /max)  as  Iwin(r,  /).  The  striations  inside  of  Iwm(r,  f  )  will  all  have 
approximately  the  same  slope  if  the  window  size  is  suff  ciently  small. 

Several  articles  have  pointed  out  the  relationship  between  the  2D-DFT  of  Iwjn(r,  /) 
and  the  slope  of  the  striations  in  /wm(r\  /)  [66,  4,  51,  18],  That  relationship  forms 
the  basis  for  the  approach  used  in  the  present  work  to  determine  the  local  striation 
angle.  The  process  of  determining  the  local  striation  angle  involves  f  ve  main  steps: 

1.  Take  a  two-dimensional  discrete  Fourier  transform  (2D-DFT)  of  7wjn(r,  /) 

2.  Eliminate  regions  of  the  2D-DFT  of  7wjn(r.  /)  associated  exclusively  with 
noise  content. 

3.  Convert  the  2D-DFT  of  7w[n(r\  /)  to  polar  coordinates. 

4.  For  several  hypothesized  striation  angles,  add  up  (integrate)  all  of  the  com¬ 
ponents  of  the  2D-DFT  of  7W in(r,  /)  corresponding  to  that  striation  angle. 

5.  The  striation  angle  that  has  the  most  “energy”  is  then  the  estimate. 
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These  steps  are  illustrated  in  Fig.  3.1  and  are  described  in  detail  in  the  follow¬ 
ing  subsections.  The  steps  are  related  to  the  Radon  Transform  and  the  Fourier- 
Slice  Theorem  [35].  If  one  were  to  skip  the  2nd  step,  then  the  steps  could  be  per¬ 
formed  with  a  Radon  Transform  using  the  Fourier-Slice  Theorem.  2  But  because 
step  2  eliminates  components  of  Iwin(r,  f)  above  particular  spatial  frequencies,  the 
present  analysis  works  directly  with  the  2D-DFT  of  IWm(r,  /). 

Interpretation  of  the  two-dimensional  discrete  Fourier  transform 

-fwin  (r, /)  has  striations  whose  slope  needs  to  be  determined  in  order  to  estimate 
the  source’s  range.  Fig.  3.1(a)  shows  an  example  /win(r,  /).  Denote  the  magnitude 
of  the  2D-DFT  of  /win(r,  /)  as  /gDF  (  At, image  •  A'y'jmage}  ■ 


hDF(kr  ,  image  5  ^/, 


image  ) 


7win(r,  f)e  j'Mfcr-.imager+fc/, image/)  df 


where  Atr, image  and  A-/amagc  do  not  refer  to  the  acoustic  wavenumbers,  but  rather 
to  the  wavenumbers  of  the  “image”  Iw in(r,f).  Hence,  Ay, image  and  At/y mage  will 
be  referred  to  as  the  image  wavenumbers.  They  are  the  horizontal  and  vertical 
axis,  respectively,  of  Fig.  3.1(b).  In  Sec.  3.3.2  the  relationship  between  image  and 
k  image  and  the  acoustic  horizontal  wavenumbers  of  the  modes  propagating  in  the 
waveguide  will  be  derived. 

In  practice  one  has  a  discrete  (sampled)  version  of  I  (r,  /),  so  the  two-dimensional 
Fourier  transform  is  implemented  as  a  two-dimensional  discrete  Fourier  transform. 
The  direction  of  the  mainlobe  originating  from  the  origin  of  bnv(kr  \m^r.,  /''/'.image) 
is  perpendicular  to  the  slope  of  the  striations  in  the  window  of  I(r,  /).  An  exam¬ 
ple  is  shown  in  Figs.  3.1(a)  and  3.1(b).  (They  do  not  appear  exactly  perpendicular 
because  of  the  different  aspect  ratios  of  the  f  gures.) 

If  one  thinks  of  Ixvm(i\  /)  as  an  image,  ignorant  of  the  fact  that  it  represents 
power  spectrum  of  an  acoustic  f  eld,  then  its  2D-DFT,  IjoviK, image,  A-/.imagc),  can 
be  interpreted  as  a  decomposition  of  /w,n(r,  /)  into  “cosine”  image  basis  functions, 
each  with  a  unique  image  wavenumber  [7,  12], 

2  The  author  thanks  an  anonymous  reviewer  for  The  Journal  of  the  Acoustical  Society  of  America 
for  pointing  this  out. 
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The  acoustician  may  gain  insight  by  noting  that  the  basis  functions  of  a  2D- 
DFT  of  an  image  look  like  two-dimensional  plane  waves  with  the  time-dependence 
removed:  exp  ( i(kxx  +  kyy)). 

Each  pixel  of  /?nr  (Ay  image  A:/, image)  represents  a  single  basis  function  with  im¬ 
age  wavenumbers  of  Ay, image  and  k  yimagc  (horizontal  and  vertical  image  wavenum¬ 
bers,  respectively,  in  Fig.  3.1(a)).  One  may  f  nd  this  easier  to  understand  if  she  or 
he  ignores  that  r  and  /  represent  range  and  frequency,  and  instead  thinks  of  them 
simply  as  labels  for  the  x  and  y  axes  of  Fig.  3.1(a). 

To  relate  image  and  A’/, image  to  striation  angles,  one  can  interpret  each  combi¬ 
nation  of  Ay,  image  and  A'/-.  iinage  (or  each  pixel  of  Fig.  3.1(b))  as  representing  a  cosine 
basis  function  at  a  particular  angle  with  a  particular  period  (striation  width). 

Application  of  the  2D-DFT  to  determine  striation  angle 

The  waveguide  invariant  as  shown  in  Eq.  (2.29)  makes  a  statement  only  about  the 
slopes  of  the  striations;  it  says  nothing  about  the  distance  between  the  striations. 
From  an  image  processing  perspective,  the  waveguide  invariant  makes  a  statement 
about  the  angle  of  the  basis  functions  comprising  Iwm(r,  /);  it  says  nothing  about 
the  period  (striation  width)  of  those  image  basis  functions.  But  as  Chuprov  points 
out  in  his  original  derivation  of  the  waveguide  invariant,  one  can  calculate  the  mini¬ 
mum  striation  width  using  only  modest  information  about  the  waveguide  [18].  This 
allows  one  to  f  Iter  out  noise  by  only  including  components  of  I(r,  /)  with  image 
wavenumbers  less  than  the  maximum  expected  due  to  the  acoustic  source  of  inter¬ 
est,  as  will  be  discussed  in  Sec.  3.3.2. 

To  determine  the  angle  of  the  striations  in  /wm(r,  /)  we  f  rst  remove  the  mean 
of  /win(r,  /)  and  then  take  its  2D-DFT  to  obtain  /2df(  Ay, image  >  A:/, image)-  Then  we 
transform  T>nF (Ay  image,  A’/amagc)  from  Cartesian  coordinates  to  polar  coordinates: 

-A2DF  (Ay  ,  image  5  ^/,  image  )=>Itd*{0,K)  (3.2) 


where 

9  =  arctan^^image/Ay^mage),  K  =  y^image  +  A;^image". 

Note  that  it  is  legitimate  to  add  Ay.jmagc  and  A: y.image  together  because  they  are  di¬ 
mensionless,  as  they  are  the  result  of  the  2D-DFT.  Flowever  throughout  most  of  the 
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present  work  Ay; image  and  Ayjmagc  are  “re-dimensionalized”  based  on  the  sampling 
used  in  the  2D-DFT,  just  as  is  typically  done  with  power  spectra  based  on  discretely 
sampled  temporal  waveforms. 

Nearest-neighbor  interpolation  is  used  to  do  the  coordinate  transform.  An  ex¬ 
ample  /2df(0i  AT),  is  shown  in  Fig.  3.1(c). 

If  one  places  bounds  on  the  possible  ranges  to  the  acoustic  source,  one  can  use 
Eq.  (2.29)  to  put  bounds  on  the  angles  that  could  have  striations  due  to  the  waveg¬ 
uide  invariant.  The  present  analysis  assumes  the  source  was  between  100  meters 
and  5,000  meters.  Typically  this  only  eliminates  a  few  degrees  (e.g.,  angles  of  3 
to  87  degrees,  instead  of  0  to  90  degrees),  but  the  striation  angle  f  nding  algorithm 
benef  ts  from  this  because  it  occasionally  incorrectly  chooses  very  steep  or  very 
shallow  angles  as  the  dominant  striation  angle. 


#min  =  arctan 


/  maximum  frequency  in  window  \ 
\  minimum  range  to  search  J 


(  minimum  frequency  in  window  \ 

C'max  =  arctan  - : - - - r — 

\  maximum  range  to  search  / 


Interpolation  isn’t  strictly  necessary  to  obtain  the  value  of  Izdf{9-,K)  at  an 
arbitrary  ( 6 ,  K ).  One  can  evaluate  the  2D-DFT  of  a  sampled  version  of  I(r,  f)  at 
arbitrary  Ay, image  and  /.'/.image  values,  analogous  to  a  discrete  time  Fourier  transform 
for  a  discrete  time  series.  However,  doing  this  is  computationally  intensive  because 
one  cannot  utilize  the  Fast  Fourier  Transform  algorithm,  and  the  present  analysis 
did  not  suggest  a  noticeable  increase  in  striation-angle  f  nding  performance  when 
doing  this  to  avoid  interpolation. 

K)  is  then  integrated  along  the  K  direction  (the  y- axis  in  Fig.  3.1(c)) 
to  add  up  all  the  components  of  the  image  with  a  particular  striation  angle  (remem¬ 
bering  that  /2df(6A  K)  will  already  have  been  spatially  fltered  to  eliminate  high 
frequency  noise). 


E[6)  =  f  I2m(e,K)dK  (3.3) 

E(9)  approximately  represents  the  amount  of  energy  (in  an  image  processing 
sense,  not  in  an  acoustic  sense)  in  huviK,  image;  &/, image)  of  striations  at  a  particular 
angle  in  /(r,  /).  A  plot  of  an  E(9)  is  shown  in  Fig.  3.1(d).  In  this  chapter,  we 
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assume  (3  =  1  and  use  E(0)  to  estimate  the  range.  However,  other  researchers 
have  modeled  (3  as  a  distribution  (see  Sec.  2.3.3)  and  used  E(6)  to  estimate  that 
distribution  [51]. 

Note  that  a  proper  change  of  variables  would  give  an  extra  factor  of  K  on  the 
right  hand  side  of  Eq.  3.3,  as  noted  in  [51].  Because  the  present  analysis  does  not 
have  that  extra  factor  of  K,  lower  values  of  K  are  weighted  more  heavily.  The 
precise  meaning  of  this  can  be  seen  by  looking  at  the  mapping  between  pixels  in 
Fig.  3.1(b)  and  Fig.  3.1(c).  This  was  done  because  it  led  to  better  estimates  of  the 
striation  angle. 

Finally,  the  angle  corresponding  to  the  maximum  value  of  E{6)  is  the  most 
dominant  angle  in  the  basis  functions  comprising  /win(r,  /).  The  angle  of  the  stri- 
ations  is  perpendicular  to  the  angle  of  the  basis  function: 

^striation  =  arg  max  E{6)  +  7r/2  (3.4) 

6 

An  huviK, image ,  /''/'.image)  obtained  from  experimental  data  will  have  noise  at 
all  values  of  Ay, image  and  Ay, image •  However,  the  region  of  /gni  ( Ay, image  *  Ay, image ) 
which  is  important  for  the  slope  estimation  (the  signal  of  interest)  lies  mostly  at 
lower  values  of  Ay:jmage  and  fc /,image  -  It  was  found  that  even  at  very  high  signal- 
to-noise  ratios,  the  noise  can  signif  cantly  deteriorate  the  slope  estimate  because 
the  maximum  value  of  Ay, image  and  Ay;image  represented  /2DF  (Ay, image:  A:/,  image)  can 
be  arbitrarily  large  depending  on  how  fnely  I(r,  f)  was  sampled  in  range  and  fre¬ 
quency.  The  effect  of  the  noise  can  be  reduced  signif  cantly  by  limiting  the  re¬ 
gion  of  integration  of  Eq.  (3.3)  to  (-Ayiimage,max  <  Ayjimage  <  Ay)image,max)  and 
(0  <  kf  image  <  k f  max)-  The  next  section  will  demonstrate  how  to  choose 

At, image, max  and  kj[mage^max. 
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Figure  3. 1 :  The  four  steps  used  to  determine  the  angle  of  the  striations  in  /win(r,  /). 
Part  (a)  shows  Iwin (r,  /),  a  window  of  I (r,  /).  The  two-dimensional  Fourier  trans¬ 
form  of  /win(r,  /)  is  then  taken  resulting  in  hmiK, image;  fc/,image)5  shown  in  part  (b) 
(using  the  bounds  described  in  Sec.  3.3.2).  hnv{kr  image.  fc/jmage)  is  then  converted 
into  polar  coordinates,  /2df($,  K),  shown  in  part  (c).  hnvW- K)  is  then  integrated 
along  K  to  produce  E (9),  shown  in  part  (d).  The  angle  corresponding  to  the  peak  of 
E{9)  is  then  the  estimated  angle  of  the  striation.  White  dotted  lines  corresponding 
to  the  estimated  striation  angle  are  then  plotted  in  part  (a)  for  visual  comparison. 
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Upper  Bounds  on  Ay,jmage  and  k /.image 


[18]  showed  that  one  can  relate  the  environmental  parameters  of  the  waveguide  to 
the  maximum  rate  at  which  I(r,  f)  can  oscillate  in  r  and  in  /.  In  this  subsection, 
we  reproduce  the  results  from  [18]  in  the  context  of  the  problem  at  hand  -  source 
range  estimation  in  a  shallow-water  waveguide.  Determining  an  upper  bound  on 
the  rate  at  which  I  (r,  f)  oscillates  in  r  and  in  /  is  equivalent  to  determining  an 
upper  bound  on  the  &v, image  and  A’/,  image  due  to  the  acoustic  source  of  interest.  Any 
components  of  72D.f(  Ay, image)  A'y'_jmage  )  above  some  Ayamage,max  and  A'/amagC_max  can 
be  regarded  as  noise  because  they  cannot  be  due  to  the  source  of  interest. 

In  order  to  exclude  as  much  noise  as  is  possible,  the  integration  in  Eq.  (3.3)  of 
-^2Df( At, image;  Aiyamagc )  will  be  bounded  by  fcr, image, max  and  Aly, image,  max-  An  example 
of  this  region  is  shown  in  Fig.  3.2.  Note  that  applying  these  bounds  on  the  is 
approximately  equivalent  to  low-pass  f  ltering  the  image. 

According  to  Eq.  (2.9),  the  acoustic  intensity  is  a  sum  of  cosines.  Each  cosine 
term  has  an  image  wavenumber  in  the  r  direction  of 


_  d  (Akmi(f)r) 

nT, image, mL  L-*rxjml  ) 

and  an  image  wavenumber  in  the  /  direction  of 

<9  (A  kmi(f)r)  rd{Akmi(f)) 


k 


/,  image,  ml 


of 


of 


(3.5) 


(3.6) 


(where  A  krnj  —  krm  kri ). 

An  upper  bound  on  Ay, image, ml  >n  /(r,  f  )  can  be  determined  (to  within  the 
approximations  used  when  deriving  the  waveguide  invariant)  by  calculating  the 
largest  possible  value  of  Akmi(f  ).  All  non-zero  values  in  /2df( Ay, image;  kf  image) 
with  Air, image  image  wavenumbers  above  the  maximum  value  of  Akm[(f)  are  likely 
due  to  noise,  and  can  be  excluded  from  the  integration  in  Eq.  (3.3).  For  all  ocean 
acoustic  waveguides,  the  horizontal  wavenumbers  are  bounded  by  [10,  Eq.  4.5.19] 


2vr/  2irf 
1 

Cmax  Cmin 


(3.7) 


where  cmax  and  cmin  are  the  respective  minimum  and  maximum  sound  speeds  that 
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kr, image  <2  71  /meters) 

Figure  3.2.  -F>df (^V, image i  kf, image)  showing  the  full  range  of  &v, image  and  kj  [mage, 
with  white  lines  illustrating  the  bounds  described  in  Sec.  3.3.2.  This  fgure  illus¬ 
trates  why  it  is  essential  to  limit  the  region  of  integration  in  Eq.  (3.3).  Even  if 
the  noise  outside  the  white  lines  is  at  a  low  level,  it  can  dominate  the  integral  in 
Eq.  (3.3)  because  of  the  large  ratio  of  outer  to  inner  areas  separated  by  the  white 
lines.  The  inner  area  is  what  is  shown  in  Fig.  3.1(b). 
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occur  in  the  environment.  Thus 


vr,  image,  max 


=  2vr  / 


(3.8) 


To  determine  and  upper  bound  on  fc/jmage.mh  we  use  Eq.  (3.6)  and  replace  / 
with  uj/(2tt),  and  note  that  is  the  reciprocal  of  the  group  speed  of  mode  m. 
The  maximum  and  minimum  group  speeds  are  approximately  bounded  by  fastest 
and  slowest  media  in  the  waveguide,  so 

kf  image, max  =  27T7’  f  J  •  (3-9) 

These  upper  bounds  are  used  to  limit  the  region  of  integration  in  Eq.  (3.3).  Us¬ 
ing  Eqs.  (3.9)  and  (3.8),  one  can  set  the  bounds  loose  enough  to  include  almost  all 
ocean  waveguides  on  Earth,  but  still  reject  much  of  the  noise  in  images  kf,  image) 

The  r  in  Eq.  (3.9)  should  be  set  to  the  maximum  range  that  one  expects  to  see  the 
source.  All  results  presented  in  this  chapter  used  Eqs.  (3.8)  and  (3.9)  with  r  =  5000 
m,  ci  =  1500  m/s,  and  C2  =  1800  m/s. 


Choosing  The  Window  Size 

Iwin(r,  /)  is  a  rectangular  window  of  I(r,  /),  bounded  by  (rmm  <  r  <  rmax)  and 
(/min  <  /  <  /max) 5  inside  of  which  the  striation  slope  will  be  estimated  using  a 
2D-DFT.  The  purpose  of  this  subsubsection  is  to  determine  how  one  goes  about 
choosing  the  window  size.  Denote  the  window  size  by 

A/  =  /max -/mm  (3.10) 

A?’  =  rmax  rm  in  (3.11) 

To  accurately  determine  the  striation  angle,  the  observation  window  of  I(r,f) 
must  be  large  enough  in  r  and  /  such  that  at  least  one  full  striation  (from  peak  to 
trough  to  peak)  is  contained  within  the  window,  in  each  direction  (r  and  /).  The 
statement  in  the  previous  sentence  can  be  quantif  ed  by  noting  that  the  frequency 
resolution  of  the  DFT  for  some  variable  x  is  A kx  =  so  if  the  DFT  of  a  signal 
is  to  distinguish  the  frequency  of  kx  from  the  zero  frequency,  then  one  needs  to 
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observe  at  least  Ax  = 

fox 

In  theory  one  could  use  the  &v, image, max  derived  Sec.  3.3.2  (originally  derived 
in  [18])  to  determine  the  minimum  Ar.  However,  only  one  term  of  the  sum  in 
Eq.  (2.11)  will  lead  to  such  a  high  value  of  Ay, image-  So  using  &v, image, max  would 
underestimate  the  minimum  Ar  that  is  required  to  estimate  the  striation  slope  in 
practical  situations.  A  better  way  to  determine  the  minimum  value  of  Ar  is  to  use 
the  value  of  &r, image, m/  averaged  over  m  and  l. 

To  do  this,  f  rst  the  approximate  values  of  A kmi  (/)  for  an  ideal  waveguide  will 
be  calculated.  Then  the  additional  approximations  required  for  non-ideal  waveg¬ 
uides  will  be  discussed.  For  an  ideal  waveguide,  the  difference  in  horizontal  wavenum¬ 
bers  of  modes  not  near  cutoff  can  be  written  as  [Eq.  (2.22)] 

AA =  km{f )  -  h(f)  «  i  Q  ^  ~  1 2)-  (3-12) 


Note  that  because  the  intensity  is  a  sum  of  cosines  and  cos(A  km,ir)  =  cos(— A  km.ir), 
only  the  absolute  value  of  A kmi  is  of  interest.  The  average  absolute  value  of  the 
wavenumber  differences  can  be  determined  by  calculating  the  average  value  of 
I  m?  —  l2\: 


( l 2  —  m 2) 


1 

Ip 


m=M  l=M 


E  El'2 

m= 1  1=1 


—  m 


M4  +  M3  -  M  '2  -  M 

Up 


(3.13) 


where  M  is  the  number  of  propagating  modes.  The  computer  algebra  system  Math- 
ematica  was  used  to  determine  the  formula  for  the  sum.  For  M  3>  1, 

(i ! 2  -  m2)  w  2  (3.14) 

3 

For  an  ideal  waveguide,  Eq.  (3.14)  can  be  inserted  into  Eq.  (3.12)  to  determine  the 
approximate  average  horizontal  wavenumber  difference,  and  thus  the  average  value 
of  the  image  wavenumber  kr  image. 

For  non-ideal  waveguides,  this  analysis  assumes  that  the  horizontal  wavenum¬ 
ber  differences  are  distributed  similarly  to  that  of  an  ideal  waveguide,  but  are 


56 


bounded  by  the  maximum  and  minimum  k  in  the  media: 


/  I  I  \  /^  —  Tflp1 

Akml(f)  =  2,f  —  -  —  W-  (3.15) 

\(-min  '-max/  1V1 

The  mean  horizontal  wavenumber  difference  is  then 

K, image, ml(f)  =  Akmt{f)  =  2t r/  (— - — ^  ^  (3.16) 

\  Cmin  ^max  /  O 

Conveniently,  this  does  not  depend  on  M. 

For  the  average  value  of  kfjmage,  we  use  a  similar  argument.  Begining  with 
Eqs.  (3.6)  and  (3.12): 


k 


/,  image,  ml 


d(A  kml(f)) 
r  df 


One  can  then  insert  M  for  an  ideal  waveguide. 

For  non-ideal  waveguides,  this  analysis  assumes  that  the  derivative  of  the  hor¬ 
izontal  wavenumbers  with  respect  to  /  are  distributed  similarly  to  that  of  an  ideal 
waveguide  (or  equivalently,  the  group  slownesses  are  distributed  similarly  to  that 
of  an  ideal  waveguide),  bounded  by  the  minimum  and  maximum  group  slownesses. 
In  that  case,  the  mean  k /, image, ml  is 


k 


/,  image,  ml 


=  27rr 


1 

Cmin 


1  'I 

3 


(3.18) 


Eqs.  (3.16)  and  (3.18)  can  be  used  to  ensure  that  the  window  will  contain  one 
full  striation  of  the  average  striation  width.  A  larger  window  size  could  be  used 
and  could  potentially  lead  to  a  more  accurate  striation  slope  estimate  because  the 
resolution  of  the  2D-DFT  is  inversely  proportional  to  the  window  size.  But  because 
the  slopes  of  the  striations  contained  in  /(r,  /)  change  with  r  and  /  according  to 
Eq.  (2.29),  the  window  size  should  not  be  too  large  or  it  will  contain  striations  with 
a  wide  range  of  slopes. 

We  now  discuss  how  to  choose  the  window  size,  given  the  trade-offs  mentioned 
in  the  previous  paragraph. 

fwin(A  /)  will  contain  striation  slopes  ranging  from  to  One  way  to 
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choose  the  window  size  would  be  to  make  the  range  of  striation  slopes  in  the  win¬ 
dow  equal  the  range  of  striation  slopes  represented  by  the  2D-DFT  bin.  According 
to  Eq.  (2.29),  a  window  of  size  A r  by  A /  centered  at  r  and  /  with  have  slopes 
ranging  from 


/  +  A//2  tQ  /  -  A//2 
r  —  Ar/2  r  +  Ar/2 


(3.19) 


A  striation’s  slope  is  perpendicular  to  the  angle  (in  r,  f  space)  of  its  cosine  basis 
function,  so  an  image  wavenumber  of  Ay, image  and  k jjmagc  represents  a  striation  with 
a  slope  of  —Ay,  image  A/,  image-  Thus,  an  image  wavenumber  frequency  bin  of  size 
AAyjmage  by  A  A'/,  image  located  at  kr^mage  and  A; /,image  represents  striations  ranging 
from 


image 


(3.20) 


To  make  the  range  of  striation  slopes  in  the  window  equal  the  range  of  striation 
slopes  represented  by  the  2D-DFT  bin,  one  could  in  principle  substitute  Ay, image  = 

and  A:/,  image  =  -fj  into  Eqs.  (3.19)  and  (3.20),  and  then  set  the  ranges  of  slopes 
equal  to  each  other.  The  solution  would  depend  on  Ay,image  and  A-jimagc,  for  which 
one  could  use  the  average  values  derived  in  this  section.  The  solution  would  also 
depend  on  r,  for  which  one  could  choose  some  value  in  the  middle  of  the  search 
range.  An  exact  analytic  solution  can  be  obtained  and  would  ensure  that  the  range 
of  slopes  represented  by  the  2D-DFT  bin  would  equal  the  range  of  slopes  in  the 
window.  Flowever,  doing  so  will  only  provide  a  relationship  between  A r  and  A/, 
not  values  for  both  quantities,  because  there  are  an  inf  nite  number  of  combinations 
of  A?’  and  A /  that  could  satisfy  the  equality. 

A  less  quantitatively  rigorous,  but  more  pragmatic  approach  is  used  in  the 
present  work  to  determine  the  value  of  A r  and  A/.  We  start  with  the  heuristic 
that  the  window  should  be  roughly  three  times  the  average  striation  width  that  we 
expect  to  see  in  each  direction.  This  heuristic  is  motived  by  a  desire  for  the  2D- 
DFT  bin  representing  the  average  image  wavenumber  expected  to  be  a  few  2D-DFT 
bin-widths  away  from  both  axes  in  Fig.  3.1(b).  Thus  we  desire: 


A  r 


3  •  2vr 

A"/,  image,  ml 


A/  = 


3  •  2vr 


image 


(3.21) 
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The  r  in  Eq.  (3.18)  should  be  set  towards  the  lower  end  of  the  ranges  over  which 
one  is  searching  for  the  source. 

In  the  present  work,  A /  was  determined  by  using  Eqs.  (3.21)  and  (3.18)  with 
ci  =  1500  m/s,  C2  =  1800  m/s,  r  =  1000  meters,  resulting  in  A /  «  81  Elz. 
Ar  was  determined  by  using  Eqs.  (3.21)  and  (3.16)  with  ci  =  1500  m/s,  C2  = 
1800  m/s,  and  /  =  525  Hz  (the  middle  of  the  frequency  range  of  the  data  presented 
in  the  next  section),  resulting  in  Ar  ~  154  meters. 

The  parameter  values  determined  in  this  section  —  Aymax,  Aymax,  A /  and  Ar 
—  are  used  to  process  the  simulated  and  both  experimental  data  sets  in  this  chapter. 
The  parameter  values  were  not  “f  ne  tuned”  for  each  data  set,  illustrating  the  ro¬ 
bustness  of  the  method  to  the  choice  of  the  parameters.  In  fact,  it  was  observed  that 
all  of  the  parameters  can  be  adjusted  by  roughly  a  factor  of  2  (increase  or  decrease) 
without  drastically  affecting  the  results  for  the  data  sets  analyzed  in  the  present 
chapter. 

Equations  (3.16),  (3.18)  and  (3.21)  provide  practical  estimates  of  the  minimum 
source  bandwidth  and  minimum  horizontal  aperture  required  for  range  estimation. 
These  estimates  are  shown  to  be  accurate  with  experimental  data  in  Sec.  3.4.2. 
However,  it  is  important  to  keep  in  mind  that  Eqs.  (3.16),  (3.18)  and  (3.21)  were  de¬ 
rived  for  range  independent  waveguides.  Under  some  circumstances  (e.g.,  a  rough 
sea  surface,  or  a  source  with  a  high  temporal  frequency  /),  effects  such  as  incoher¬ 
ent  scattering  will  smear  out  the  high  image-wavenumber  striations,  and  a  window 
larger  than  that  predicted  by  Eq.  (3.21)  will  be  necessary  to  estimate  a  striation’s 
slope.  (See  Sec.  7.3  for  specif  c  suggestions  for  future  research) 

3.4  Results 

3.4.1  2D-DFT  technique  applied  to  simulated  data 

In  this  subsection,  the  striation  angle  f  nding  technique  based  on  the  2D-DFT  is 
tested  on  the  simulated  acoustic  intensity  shown  in  Fig.  3.3,  which  is  from  the 
same  Pekeris  waveguide  used  to  produce  Fig.  2.3.  After  f  nding  the  striation  angle, 
Eq.  (3.1)  is  used  to  estimated  the  range  to  the  acoustic  source,  assuming  (3  =  1. 

First,  I(r,  f )  is  divided  into  several  windowed  segments,  IWm{r,  /),  spread  on 
a  grid  throughout  the  (r,  /)  plane.  The  striation  angle  is  then  estimated  in  each 
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Figure  3.3:  Acoustic  intensity  (dB,  arbitrary  reference)  in  a  Pekeris  waveguide 
plotted  versus  range  and  frequency,  I(r,f),  for  a  receiver  depth  of  20  meters  and  a 
source  depth  of  40  meters. 
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/wm  (r,  f)  using  the  technique  described  in  the  previous  section. 

The  resulting  slope  estimates  are  shown  in  Fig.  3.4(a).  Each  slope  estimate  is 
associated  with  a  range  estimate,  which  are  shown  in  Fig.  3.4(b). 

In  order  to  use  all  of  the  data,  one  needs  to  ensure  that  every  pixel  in  Fig.  3.3  is 
contained  inside  at  least  one  window.  In  fact,  one  can  let  the  windows  overlap.  In 
the  present  analysis,  the  windows  overlapped  by  about  50%  in  both  the  r  and  the  / 
directions. 

One  can  average  each  column  of  Fig.  3.4(b)  to  obtain  more  accurate  range 
estimates.  The  result  of  such  averaging,  plotted  versus  the  true  range  to  the  source, 
is  shown  in  Fig.  3.4(c). 

The  estimates  are  biased  by  about  10%.  This  is  most  likely  do  to  the  fact  that  the 
derivation  of  the  waveguide  invariant  does  not  take  into  account  the  £  cylindrical 
spreading,  which  will  cause  the  actual  slopes  to  be  steeper  than  that  predicted  by 
the  waveguide  invariant.  However  it  could  also  be  due  to  the  other  approximations 
made  when  the  waveguide  invariant,  such  as  the  dependence  of  the  mode  shapes  on 
frequency.  Because  this  error  is  less  than  the  expected  range  estimate  accuracy  for 
experimental  data,  it  will  not  be  addressed  further  in  the  present  analysis. 

3.4.2  2D-DFT  technique  applied  to  experimental  data 

The  experimental  data  presented  in  this  chapter  was  collected  during  GLINT08,  an 
experiment  performed  during  the  summer  of  2008  near  Pianosa  Island,  Italy.  Two 
sound  speed  prof  les  measured  about  an  hour  before  and  an  hour  after  the  acoustic 
data  were  collected  are  shown  in  Fig.  3.5(a),  but  note  that  this  information  was  not 
used  by  the  range  estimation  algorithm. 

An  acoustic  source  was  lowered  40  m  below  the  ocean  surface  from  the  R/V 
Leonardo,  which  was  using  dynamic  positioning  to  keep  its  position  as  f  xed  as 
possible.  Due  to  a  malfunctioning  GPS  unit,  the  acoustic  source’s  position  had  an 
uncertainty  of  100  m.  The  signal  projected  from  the  source  was  pseudo-random 
white  noise  with  a  with  an  approximately  f  at  spectrum  from  300  to  750  Hz  and  a 
frequency-integrated  source  level  of  150  dB  re:  1  //Pa  at  1  m.  Due  to  equipment 
limitations,  the  signal  was  a  single  65000  point  realization  of  white  noise  that  was 
repeated  back-to-back  every  5.33  seconds. 

An  acoustic  receiver  was  towed  directly  toward  the  acoustic  source  from  a  range 
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(a)  Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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(c)  Estimated  Ranges  vs.  True  Ranges 
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Figure  3.4:  Range  estimates  for  a  simulated  Pekeris  waveguide,  (a)  Same  as 
Fig.  3.3,  but  with  the  estimated  slopes  superimposed  as  white  solid  lines,  (b)  Range 
estimates  based  on  each  slope  in  part  (a)  plotted  versus  the  true  range,  (c)  Range 
estimates  versus  true  range,  obtained  by  averaging  each  column  of  part  (b).  The 
bias  is  discussed  Sec.  3.4.1. 
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(a)  Sound  Speed 
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Figure  3.5:  (a)  Sound  speed  prof  les  taken  before  and  after  the  acoustic  data  col¬ 
lection.  (b)  Sound  speed  at  the  hydrophone’s  depth  plotted  versus  range  along  the 
track,  (c)  Water  column  depth  (solid  line)  and  hydrophone  depth  (dashed  line) 
plotted  versus  range  from  the  source  along  the  hydrophone’s  path. 
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of  2,200  meters  to  about  500  meters  at  a  speed  of  1.5  m/s  and  at  a  depth  of  30 
meters.  This  segment  of  the  experiment  is  referred  to  as  the  incoming  segment. 
The  acoustic  receiver  was  then  lowered  to  50  meters,  and  towed  away  from  the 
source  at  1.5  m/s  back  out  to  a  range  of  2,200  meters.  This  segment  is  referred  to 
as  the  outgoing  segment.  The  acoustic  receiver  was  moving  continuously,  so  the  all 
of  the  presented  data  was  collected  in  less  than  1  hour.  There  was  a  small  Doppler 
shift  of  about  1%  due  to  receiver  motion.  Such  a  shift  may  cause  range  estimate 
errors  of  about  1%,  which  is  negligible  compared  to  the  overall  expected  accuracy 
and  thus  will  be  ignored.  The  acoustic  receiver  location  had  an  uncertainty  of  50 
meters. 

The  acoustic  data  were  sampled  at  a  frequency  of  4  kHz.  The  window  length 
used  to  estimate  the  spectrum  of  the  received  signal  had  to  be  21,334  points  long  so 
that  it  corresponded  to  the  5.33  second  repetition  rate  of  the  signal.  Under  other  cir¬ 
cumstances,  one  would  be  free  to  choose  other  window  lengths  that  do  not  strongly 
depend  upon  the  signal  characteristics.  Because  the  acoustic  receiver  was  moving 
at  an  approximately  1.5  m/s,  each  spectrum  (each  column  of  Fig.  3.6)  represents 
~  5.33  s  x  1.5  m/s  =  8  meters  in  distance  that  the  acoustic  receiver  traveled.  The 
Blackman-Tukey  method  of  spectrum  estimation  was  used  to  estimate  each  spec¬ 
trum  (each  column  of  Fig.  3.6).  Because  the  acoustic  receiver  was  towed  at  an 
approximately  constant  rate,  the  spectrogram  of  the  recorded  time  series  is  /(r,  /). 

Once  the  spectrogram,  I(r,  /),  was  calculated,  the  processing  method  was  ex¬ 
actly  the  same  as  that  used  for  the  simulated  data.  The  same  parameter  values  were 
used. 

Incoming  segment 

The  water  column  depth,  acoustic  receiver’s  depth  and  sound  speed  along  the  hy¬ 
drophone’s  path  are  plotted  in  Fig.  3.5(b)  and  (c).  The  measured  acoustic  f  eld  from 
0  Hz  to  2  kHz  is  shown  in  Fig.  3.6.  The  results  of  the  range  estimation  algorithm 
are  shown  in  the  same  format  as  the  simulated  results,  in  Figs.  3.7(b)  and  3.7(c). 

There  is  extremely  good  qualitative  agreement  between  the  angle  determined 
by  the  striation  angle  f  nding  algorithm  and  the  striation  angle  as  it  appears  to  a 
human  observer,  as  can  be  seen  in  Fig.  3.7(a).  At  source  ranges  larger  than  1000 
meters,  the  estimated  range  tends  to  be  less  than  the  true  range.  This  could  be  the 
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result  of  any  number  of  effects,  including:  3-d  propagation  effects,  range  inhomo¬ 
geneities,  the  temporal  stationarity  of  the  SSP  while  the  f  eld  was  being  measured, 
the  approximations  used  when  deriving  the  waveguide  invariant.  Sec.  3.4.3  shows 
a  simulated  spectrogram  for  comparison. 
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Figure  3.6:  The  acoustic  intensity  (dB,  arbitrary  reference)  as  measured  by  the 
acoustic  receiver  for  the  incoming  segment  of  the  experiment.  The  striations  are 
clearly  visible  in  the  frequency  range  of  the  source,  from  350  to  700  Flz. 
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Outgoing  Segment 

The  watercolumn  depth  and  sound  speed  along  the  hydrophone’s  path  were  similar 
to  the  incoming  segment  (see  Fig.  3.5).  The  measured  acoustic  feld,  from  0  to 
2  kHz  is  shown  in  Fig.  3.8.  There  is  a  loud  interfering  broadband  acoustic  source 
during  the  portion  of  the  spectrogram  corresponding  to  a  range  of  1000  to  2000 
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(a)  Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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Figure  3.7:  Range  estimates  for  the  incoming  segment  of  the  experiment,  (a)  A 
zoom-in  of  Fig.  3.6  on  the  frequencies  of  interest,  with  the  estimated  slopes  super¬ 
imposed  as  solid  white  lines,  (b)  Range  estimates  based  on  each  estimated  slope 
in  part  (a),  plotted  versus  the  true  range,  (c)  Range  estimates  versus  true  ranges, 
obtained  by  averaging  each  column  of  part  (b). 
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meters.  As  one  can  see  in  Fig.  3.9,  this  does  effect  the  estimates,  but  not  as  much 
as  one  might  expect.  The  estimates  are  only  adversely  affected  at  ranges  of  1400  to 
1 800  meters,  when  the  striations  in  the  spectrogram  from  the  interferer  are  nearly 
parallel  to  the  striations  that  are  expected  from  the  experimental  acoustic  source. 
If  one  did  not  know  a  priori  that  the  acoustic  receiver  was  moving  away  from  the 
acoustic  source  of  interest,  then  the  estimates  may  have  been  affected  more  ad¬ 
versely  because  one  would  have  had  to  search  over  the  full  180  degrees.  The  range 
estimates  are  accurate  within  a  few  hundred  meters,  excluding  the  estimates  when 
the  true  range  was  1400  to  1800  meters. 


Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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Figure  3.8:  The  acoustic  intensity  (dB,  arbitrary  reference)  as  measured  by  the 
acoustic  receiver  along  its  track  away  from  the  source.  Note  the  striations  in  the 
frequencies  of  the  acoustic  source  (350  to  700  Hz)  and  the  interfering  source  present 
during  times  corresponding  to  ranges  from  1000  to  2000  meters. 
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(a)  Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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range  (m) 


Figure  3.9:  Range  estimates  for  the  outgoing  segment  of  the  experiment,  (a)  A 
zoom-in  of  Fig.  3.8  on  the  frequencies  of  interest,  with  the  estimated  slopes  super¬ 
imposed  as  solid  white  lines,  (b)  Range  estimates  based  on  each  estimated  slope 
in  part  (a),  plotted  versus  the  true  range,  (c)  Range  estimates  versus  true  range, 
obtained  by  averaging  each  column  of  part  (b).  Large  errors  occur  when  the  inter¬ 
fering  source  was  present. 
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So  far,  the  present  chapter  has  dealt  exclusively  with  the  signal  recorded  from 
a  single  hydrophone.  However,  data  was  collected  from  a  32-element  horizontal 
array  of  hydrophones.  Chapter  5  discusses  how  to  process  the  data  from  all  the 
hydrophones  in  order  to  reduce  the  noise  in  Fig.  3.8.  The  result  of  such  f  ltering 
can  be  seen  in  Fig.  5.4(b).  One  can  then  apply  the  passive  ranging  algorithm  to 
the  array-f  ltered  data  in  Fig.  5.4(b);  the  result  is  shown  in  Fig.  3.10.  Note  that 
Figs.  3.8  and  3. 10  are  from  the  experimental  data  set;  the  difference  is  that  Fig.  3.8 
used  data  from  a  single  hydrophone  of  the  array  whereas  Fig.  3.10  used  data  all  of 
the  hydrophones  of  the  array. 
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(a)  Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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Figure  3.10:  Range  estimates  for  the  outgoing  segment  of  the  experiment  using 
array-f  ltered  data.  This  f  gure  is  similar  to  Fig.  3.9,  but  data  from  the  32-element 
array  has  been  used  to  remove  the  noise  from  sources  besides  the  source  of  interest, 
(a)  A  zoom-in  of  Fig.  5.4  on  the  frequencies  of  interest,  with  the  estimated  slopes 
superimposed  as  solid  white  lines,  (b)  Range  estimates  based  on  each  estimated 
slope  in  part  (a),  plotted  versus  the  true  range,  (c)  Range  estimates  versus  true 
range,  obtained  by  averaging  each  column  of  part  (b). 
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3.4.3  Validity  of  (3  =  1  assumption 

Section  3.2  discussed  that  although  usually  (3  «  1  in  shallow  water  environments, 
that  may  not  always  be  the  case  (see  Sec.  2.3.3  or  [51]).  Determining  when  one  can 
assume  that  (3  «  1  is  still  an  active  area  of  research,  but  some  guidelines  are  given 
in  [51,  49],  and  other  papers  on  the  waveguide  invariant.  Chapter  4  presents  some 
original  research  on  this  topic. 

All  of  the  analysis  in  the  present  chapter  assumed  (3  =  1,  so  it  is  worthwhile 
to  simulate  the  acoustic  f  eld  in  an  environment  similar  to  the  environment  where 
the  experimental  data  was  collected,  in  order  to  determine  if  (3  differs  signif  cantly 
from  the  assumed  value.  (We  expect  it  not  to,  otherwise  we  wouldn’t  have  been 
able  to  accurately  estimate  the  range  to  the  source  in  the  previous  section.) 

To  do  this,  the  normal  mode  program  Kraken  [48]  was  used.  The  sound  speed 
prof  le  used  was  that  collected  during  the  experiment  [shown  in  Fig.  3.5(a)],  The 
source  and  receiver  geometry  were  the  same  as  the  incoming  portion  of  the  exper¬ 
iment.  Because  the  bottom  bathymetry  and  the  sound  speed  prof  le  changed  very 
little  with  range,  the  environment  was  modeled  as  being  range  independent.  The 
bottom  properties  were  unknown  so  typical  values  for  a  bottom  half-space  were 
used  (Cbottom  =  1650 m/s,  bottom  =  1.5g/cm3,  a  =  0.5dB/A). 

Figure  3.11  shows  the  simulated  spectrogram.  The  white  lines  have  slopes 
corresponding  to  (3  =  1.  It  can  be  seen  that  there  are  some  striations  with  slopes 
that  are  slightly  steeper  than  the  (3  =  1  lines,  but  almost  no  striations  have  slopes 
that  are  shallower  than  the  (3  =  1  lines.  This  could  explain  why  the  ranges  tended 
to  be  underestimated  in  some  parts  of  the  incoming  segment  of  the  experiment 
[Fig.  3.7(c)],  and  suggests  that  (3  may  have  a  value  slightly  larger  than  one  for  this 
particular  environment  and  source-receiver  geometry. 

The  fact  that  the  range  estimates  from  the  experimental  data  in  Figs.  3.7(c) 
and  3.10(c)  were  between  about  75%  to  100%  of  the  true  range  implies  that  (3 
had  a  value  between  about  1  and  4/3  for  the  environment  where  the  experimental 
data  was  collected.  Later  in  this  thesis,  a  method  is  developed  for  calculating  the 
value  of  (3.  Section  4.5.4  shows  that  (3  changes  with  horizontal  wavenumber,  and 
that  it  has  a  value  of  slightly  less  than  one  for  this  particular  environment  (for 
the  SRBR  modes).  Thus  there  is  a  discrepancy  between  the  predicted  value  of  (3 
being  slightly  less  than  one,  and  the  measured  value  of  f3  being  slightly  greater  than 
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one  ([49]  and  [32]  also  have  experimental  examples  of  /3  being  slightly  greater 
than  one).  This  discrepancy  might  be  due  to  the  non-SRBR  modes  which  have 
different  fj  values,  the  1  / \fr  spreading,  or  attenuation  (imaginary  components  of 
the  horizontal  wavenumbers).  [51]  provides  some  explanation  for  this  discrepancy, 
but  there  is  certainly  room  for  further  research. 


Acoustic  Intensity  (dB,  Range  vs.  Frequency) 
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Figure  3.11:  A  simulated  spectrogram  for  the  environment  used  to  collect  the  ex¬ 
perimental  data.  The  white  lines  correspond  to  /3  =  1.  Most  of  the  striations 
correspond  to  /3  «  1,  with  a  few  exceptions. 


3.5  Chapter  summary 

A  processing  scheme  based  on  the  waveguide  invariant  and  the  2D-DFT  of  I(r,  /) 
was  used  to  estimate  an  acoustic  source’s  range  using  simulated  data  and  two  sets 
of  experimental  data.  The  processing  techniques  used  did  not  require  human  inter- 
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pretation  of  any  images  in  order  to  obtain  the  range  estimate,  making  the  techniques 
suitable  for  implementation  on  an  autonomous  platform. 

A  relationship  between  the  average  image-wavenumbers  in  the  2D-DFT  of 
/ (r,  /)  and  the  acoustic  waveguide  parameters  was  used  to  determine  the  mini¬ 
mum  observation  window  size  of  I(r,  f)  required  for  range  estimation.  A  similar 
relationship  was  used  to  reject  noise  in  I(r,  /). 

The  same  set  of  signal  processing  parameter  values  (  maximum  image  wavenum¬ 
ber,  and  window  size  in  range  and  frequency)  was  used  for  both  simulated  and  ex¬ 
perimental  data,  showing  that  the  signal  processing  parameters’  values  do  not  need 
to  be  f  ne-tuned  for  each  data  set.  The  range  estimates  were  based  on  the  assump¬ 
tion  (approximation)  that  (5  is  a  single  scalar  value,  and  is  equal  to  one.  The  range 
estimates  were  accurate  to  within  about  25%  despite  using  only  a  minimal  amount 
of  a  priori  knowledge  about  the  environment.  The  technique  appeared  to  be  robust, 
but  more  experimental  data  would  be  needed  to  determine  the  robustness  of  these 
algorithms  in  other  environmental  conditions. 
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Chapter  4 

A  WKB  modal  approach  to 
calculating  the  waveguide 
invariant 


The  previous  chapter  illustrated  that  the  waveguide  invariant  is  a  powerful  con¬ 
cept  for  explaining  range-frequency  striations  in  experimental  data.  It  is  somewhat 
surprising  that  the  strong  non-uniformities  present  in  experimental  SSP  in  the  pre¬ 
vious  chapter  did  not  have  a  large  effect  on  the  value  of  f3  observed  in  I(r,  /).  The 
present  chapter  develops  a  method  for  calculating  the  waveguide  invariant  that  il¬ 
lustrates  the  relationship  between  the  SSP  and  /3,  which  provides  some  insight  into 
why  /3  «  1  in  most  shallow-water  waveguides.1 

In  this  chapter,  m  and  l  refer  to  the  mode  numbers  just  as  they  did  in  previous 
chapters.  This  chapter  also  uses  n  to  refer  to  the  mode  number,  but  n  is  used  only 
when  the  mode  number  comes  from  the  WKB  approximation.  The  reason  for  this 
distinction  will  become  clear  later. 


'This  chapter  is  based  on  “A  modal  Wentzel-Kramers-Brillouin  approach  to  calculating  the 
waveguide  invariant  for  non-ideal  waveguides”  by  Kevin  L.  Cockrell  and  Henrik  Schmidt  [Under 
review  for  publication  in  The  Journal  of  the  Acoustical  Society  of  America] 
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4.1  Introduction 


In  Chuprov’s  original  paper  on  the  waveguide  invariant,  he  def  nes  j3  in  two  differ¬ 
ent  ways.  The  f  rst  way  uses  normal  modes  to  describe  the  acoustic  f  eld,  and  the 
second  way  uses  ray  theory  to  describe  the  acoustic  f  eld  [18].  In  general,  normal 
modes  is  accurate  under  a  wider  range  of  circumstances  than  ray  theory,  however, 
normal  modes  tends  to  be  less  amenable  to  non-numerical  analysis  than  ray  theory. 

Previous  research  on  the  waveguide  invariant  using  ray  theory  has  used  both 
numerical  techniques  and  analytic  techniques  [18,  24,  22],  with  the  analytic  tech¬ 
niques  providing  some  insight  into  how  the  SSP  affects  j3. 

Previous  research  on  the  waveguide  invariant  using  normal  modes  has  mostly 
used  numerical  techniques  to  calculate  (3  \  18,51,21,  52],  which  has  provided  many 
useful  insights  but  is  not  conducive  for  understanding  how  the  SSP  affects  (3. 

This  chapter  investigates  the  waveguide  invariant  using  a  the  WKB  approxima¬ 
tion  to  normal  modes,  which  can  be  thought  of  as  “half  way”  between  full-f  eld 
normal  modes  and  geometric  ray  theory.  The  modal  WKB  description  is  accu¬ 
rate  under  a  wider  range  of  circumstances  than  geometric  ray  theory  but  is  more 
amenable  to  (non-numerical)  analysis  than  full-f  eld  normal  modes. 


/3  can  be  def  ned  directly  in  terms  of  the  normal  modes’  horizontal  wavenum¬ 
bers,  as  is  shown  in  Eq.  (2.15).  To  make  the  dependence  of  (3  on  the  modes  numbers 
m  and  l  and  the  frequency  oo  more  explicit,  we  write  Eq.  (2.15)  as: 


/ 3{m ,  l ,  to ) 


1 

UJ 


kr(m ,  ui)  —  kr(l ,  to) 


d(kr(m,u)—kr(l,w)) 

():jJ 


(4.1) 


where  kr(i ,  u)  is  the  horizontal  wavenumber  of  mode  i  at  frequency  c o. 

In  Sec.  2.3.2  it  was  shown  analytically  that  (3  «  1  for  ideal  waveguides.  Ref¬ 
erences  [18]  and  [11]  analytically  show  that  f3  ~  —3  for  surface-trapped  modes  in 
waveguides  with  n2- linear  sound  speed  prof  les  (SSP).  As  discussed  in  Sec.  2.3,  ex¬ 
perimental  data  and  numerical  simulations  have  shown  that  many  realistic  shallow- 
water  waveguides  can  be  well  approximated  by  an  ideal  waveguide  for  the  purpose 
of  calculating  (3  (i.e.,  (3  «  1  for  many  shallow-water  waveguides).  However,  the 
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analysis  in  Sec.  2.3  (originally  from  [18],  [11,  Sec.  6.7.2],  and  [34,  Sec.  2.4.6]) 
does  not  lend  itself  to  understanding  why  B  is  not  strongly  affected  by  the  non¬ 
uniformities  present  in  a  realistic  SSPs,  as  opposed  to  the  pressure  f  eld  itself  which 
is  strongly  affected  (two  examples  are  shown  in  Fig.  2.4).  In  fact,  the  derivation  in 
Sec.  2.3.2  showing  that  /3  «  1  for  an  ideal  waveguide  only  considers  modes  far 
from  cutoff  (i.e.,  low-order  modes  at  high  frequency).  Paradoxically,  those  are  pre¬ 
cisely  the  modes  whose  kr  values  will  be  most  affected  by  the  non-uniformities  in 
the  SSP  2;  those  modes  will  “tuck”  into  the  local  sound  speed  minimums  and  will 
not  behave  like  modes  in  an  ideal  waveguide. 

This  chapter  develops  a  method  for  calculating  B  based  on  the  WKB  approxi¬ 
mation,  which  is  an  approximate  method  for  solving  for  the  horizontal  wavenum¬ 
bers  in  a  waveguide.  The  relationship  between  the  WKB  approximation  and  ft  is 
brief  y  mentioned  in  the  original  waveguide  invariant  paper  [18].  Later,  Brown  et 
al.  used  the  WKB  approximation  to  show  the  relationship  between  between  f3  and 
the  ray  stability  parameter  a  [14,  6].  The  present  chapter  investigates  how  non¬ 
uniformities  in  a  SSP  will  affect  (5.  The  resulting  method  is  then  applied  to  a  few 
canonical  waveguides  and  compared  with  the  values  of  B  calculated  numerically 
from  the  normal  mode  program  Kraken  [48], 

The  main  results  are: 

•  If  the  WKB  approximation  is  used,  the  derivatives  required  to  calculate  B 
can  be  performed  implicitly,  circumventing  the  need  to  obtain  explicit  so¬ 
lutions  for  kr{m,u)  and  allowing  for  the  inclusion  of  a  fuid  bottom  halfs¬ 
pace.  The  resulting  estimate  of  B  is  then  an  explicit  function  of  the  horizontal 
wavenumber  instead  of  being  an  explicit  function  of  the  mode  number. 

•  When  the  bottom  halfspace  is  modeled  as  a  vacuum,  the  value  of  B  for  all 
horizontal  wavenumbers  and  all  frequencies  is  a  function  of  a  single  param¬ 
eter:  the  phase  slowness  Sp. 

•  The  value  of  B  can  be  directly  related  to  two  one-dimensional  curves:  the 
depth-integrated  vertical  wavenumber  as  a  function  of  phase  slowness  and 
the  seaf  oor  ref  ection  coeff  cient  as  a  function  of  phase  slowness.  These 

2Chuprov  mentions  this  fact  later  in  his  paper  when  discussing  the  waveguide  invariant  in  deep¬ 
water  waveguides. 
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curves  show  that  /3  for  surface-ref  ected  bottom-ref  ected  (SRBR)  modes  is 
not  strongly  affected  by  the  details  of  the  SSR,  but  that  (3  for  non-SRBR 
modes  is  strongly  affected  by  the  details  of  the  SSR 


•  A  qualitative  argument  suggests  that  (3  k,  1  for  SRBR  modes  in  waveguides 
where  the  seaf  oor  is  modeled  as  a  vacuum,  regardless  of  the  SSR  (3  can  take 
on  a  wide  range  of  values  for  non-SRBR  modes. 


4.2  Calculating  (3  using  the  WKB  approximation 


To  calculate  f3,  the  WKB  approximation  is  used  to  determine  the  modal  horizontal 
wavenumbers.  Details  of  the  WKB  approximation  as  it  applied  to  calculating  the 
modal  horizontal  wavenumbers  can  be  found  in  [11,  Sec.  6.7],  [33,  Sec.  2.5],  [5], 
and  [60,  Sec.  2.9]  (See  Appendix  B  for  even  more  references).  All  results  in  this 
chapter  apply  only  to  propagating  modes  with  no  more  than  2  turning  points  (depths 
where  kz(z,kr)  =  0.  see  Sec.  2.2.2),  and  to  waveguides  with  bottoms  that  are 
either  lossless  homogeneous  f  uid  halfspaces  or  vacuums.  Because  the  quantity 
of  interest  in  the  present  analysis  is  the  acoustic  intensity,  we  note  that  the  WKB 
approximation  is  usually  more  accurate  at  calculating  the  acoustic  intensity  than 
it  is  at  calculating  the  complex  acoustic  f  eld.  As  discussed  in  [8,  Sec.  48.5],  this 
is  because  the  error  of  the  WKB-approximated  acoustic  intensity  depends  on  the 
errors  of  the  differences  of  the  horizontal  wavenumbers,  as  opposed  to  the  error  of 
the  WKB-approximated  complex  acoustic  f  eld  which  depends  on  the  errors  of  the 
value  of  the  individual  horizontal  wavenumbers. 

Under  the  WKB  approximation,  the  horizontal  wavenumbers  are  calculated  by 
solving  for  Ay  in  the  equation  [5] 

n  =  (j)(kr,u )  +  A0dn(ky,  oj)  +  A</)up(Ay,  u)  +  1  (4.2) 
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where  n  is  the  mode  number  (1,  2,  3,  ...), 


1  f z 2 

< j)(kr,uj )  =  —  /  fcz(z,  kr)  dz  (4.3) 


where  zi  and  z-i  are  either  the  sea  surface  (z  =  0)  and  sea  f  oor  (z  =  d)  respectively, 
or  a  z  at  which  kz(z)  =  0  (a  turning  point).  And  with: 

if  kr  <  k( 0)  (i.e.,  mode  ref  ects  from  surface) 

if  kr  >  k( 0)  (i.e.,  mode  has  an  upper  turning  point), 

(4.5) 


if  kr  <  k(d )  (i.e,  mode  ref  ects  from  bottom) 

if  kr  >  k(d)  (i.e.,  mode  has  a  lower  turning  point), 

(4.6) 

where  4>b(kr,  u)  is  jz  times  the  phase  angle  of  the  bottom  half-space  refection 
coeffcient.  (—1/2  <  cf>b(kr,u))  <  0)  for  a  fuid  bottom  and  <f>b(kr,uj)  =  —1/2  for 
a  vacuum  bottom.  Figure  4.1  discusses  how  to  interpret  Eq.  (4.3)  (See  Fig.  2.1  for 
a  def  nition  of  SRBR  modes). 


^4>dn(kr,Uj')  — 


4>b 
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Figure  4.1:  Illustration  of  the  WKB  approximation  applied  to  calculating  the 
horizontal  wavenumbers  kr.  The  shaded  area  is  related  to  (but  not  equal  to) 
f  kz (z,  kr )  dz,  which  is  used  to  determine  the  values  of  kr  that  correspond  to  propa¬ 
gating  modes.  Subf  gure  (a)  demonstrates  that  when  kr  is  greater  than  the  minimum 
k(z)  in  the  water  column  (which  corresponds  to  non-SRBR  modes),  the  shape  of 
f  kz(z,  kr)  dz  as  a  function  of  kr  will  depend  strongly  on  the  detailed  shape  of  the 
the  SSP  —  especially  the  part  of  the  SSP  near  k{z)  =  Ay.  Subf  gure  (b)  demon¬ 
strates  that  when  kr  is  less  than  the  minimum  k[z)  in  the  water  column  (which 
corresponds  to  SRBR  modes),  the  shape  of  f  kz(z ,  Ay)  dz  as  a  function  of  kr  does 
not  depend  strongly  on  the  details  of  the  SSP  because  the  integral  tends  to  average 
out  any  “roughness”  in  the  SSP.  A  similar  argument  can  be  made  regarding  the 
dependence  of  J  kz(z )  dz  on  co. 
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Eq.  (2.15)  can  be  manipulated  to  facilitate  a  geometric  interpretation  of  f3  in 
terms  of  the  WKB  approximation.  The  numerator  and  denominator  of  Eq.  (2.15) 
can  be  divided  by  an  arbitrary  value  without  affecting  the  value  of  the  fraction. 
Specif  cally,  one  can  divide  the  numerator  and  denominator  by  the  difference  in  the 
mode  numbers  ( m  —  l )  so  that 


1  (kr(m,u>)  —  kr(l,uj))  /(m  —  l) 

UJ  d(kr(m,u)-kr(l,u>))  j  ^  ^ 


(4.7) 


A  line  connecting  the  points  ( kr(m ,  ui),  m)  and  (kr(l,  u),  l),  as  shown  in  Fig.  4.2 
with  m  =  2  and  l  =  4,  has  a  slope  of  (m  —  l)/(kr(m ,  u)  —  kr(l,  u)).  This  allows 
one  to  interpret  the  numerator  of  Eq.  (4.7)  as  the  reciprocal  of  that  slope. 


(a)  (b) 


Figure  4.2:  (a)  A  plot  of  n  versus  kr  for  an  ideal  waveguide  using  the  WKB  approx¬ 
imation  (/  =  100  Flz,  d  =  50  m,  c  =  1500  m/s).  The  two  points  plotted  correspond 
to  modes  2  and  4.  The  reciprocal  of  the  slope  of  the  line  connecting  the  two  points 
can  be  used  to  facilitate  interpretation  of  the  waveguide  invariant  /3,  as  described  in 
Sec.  4.2.  (b)  A  plot  of  n  versus  kr  for  a  waveguide  with  an  n2 -linear  sound  speed 
prof  le  and  pressure-release  boundary  conditions  (/  =  100  Flz,  d  =  50  m,  c(0)  = 
1400  m/s,  c(d)  =  1600  m/s).  The  discontinuity  at  kr  =  u/c(d )  is  due  to  the  WKB 
approximation  (//yn  as  def  ned  in  Eq.  (4.6)). 

/3(m,  l,u> )  depends  on  how  ( kr(m,u )  —  kr(l,oj ))  depends  on  frequency.  One 
can  visualize  that  dependence  by  plotting  n  as  a  function  of  kr  and  of  u,  as  shown 
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for  an  ideal  waveguide  in  Fig.  4.3. 


n  vs.  k  and  co 

r 


k  (1/m) 


Figure  4.3:  A  contour  plot  of  n  versus  kr  and  c o  for  an  ideal  waveguide  (d  =  50  m, 
c  =  1500  m/s).  The  solid  black  lines  are  contours  at  integer  values  of  n,  which 
in  the  WKB  approximation  correspond  to  the  normal  modes’  horizontal  wavenum¬ 
bers. 

/3(m,  l,u)  can  then  be  interpreted  as:  -1  divided  by  the  slope  of  the  line  con¬ 
necting  the  points  (/cr,cu,m)  and  ( kr,u,l ),  divided  by  how  the  reciprocal  of  that 
slope  changes  with  frequency.  Or, 


1/slope 


<9(l/slope)  \ 

duj  )  i 
/  m,l 


(4.8) 


where  the  m,  l  subscript  indicates  that  the  m  and  l  values  at  which  the  points  are 
located  must  stay  constant  as  u  changes. 

The  discussion  in  the  caption  of  Fig.  4.1  suggests  that  (kr(m,u)  —  kr(l,u )) 
and  its  dependence  on  frequency  does  not  strongly  depend  on  the  details  c(z)  when 
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| kr(m,  u)  —  kr(l,  w)|  <C  fcmin  and  kr(m,  u)  is  not  too  close  to  km ;n,  in  some  quali¬ 
tative  sense.  So  then  one  would  expect  that  /?  would  also  not  depend  on  the  details 
of  the  SSP  under  those  circumstances  in  waveguides  with  a  pressure-release  bot¬ 
toms.  This  is  mentioned  by  Chuprov  in  [18,  pp.  104], 


4.2.1  Approximating  a  finite  difference  with  a  continuous  derivative 

The  acoustic  f  eld  in  a  waveguide  is  typically  dominated  by  a  few  groups  of  modes 
with  close  mode  numbers,  for  reasons  discussed  in  Sec.  2.2.3.  Because  of  that, 

/3  is  often  defned  only  in  terms  of  adjacent  modes  (see  Sec.  2.3.3).  Because  the 
mode  number  n  is  a  monotonic  function  of  the  horizontal  wavenumber  kr,  a  group 
of  adjacent  modes  corresponds  to  a  particular  range  of  kr  values.  Here,  we  use 
that  relation  to  approximate  a  f  nite  difference  in  mode  order  with  a  continuous 
derivative  in  n. 

Ifthe  curve  n(kr,u)  is  approximately  linear  in  kr  between  the  points  (n,  kr(n  ,  u)) 
and  (n  +  An,  kr  (n  +  An  ,  w)),  then  the  f  nite  difference  between  horizontal  wavenum¬ 
bers  in  Eq.  (4.7)  can  be  approximated  by  a  continuous  partial  derivative: 

f  (kr(n  +  An,w)  -  kr(n,u))\  _  dkr(n,uj) 

A™o  V  An  )  ~  dn  '  1  '  ’ 

This  is  equivalent  to  the  idea  of  expanding  the  group  and  phase  slownesses  around 
an  average  value  for  modes  with  adjacent  kr  values,  as  discussed  in  Sec.  2.3.3,  [18], 
and  [11,  Sec.  6.7.2],  Most  of  the  present  analysis  applies  only  to  mode  pairs  with 
adjacent  kr  values,  but  Sec.  4.6  discusses  how  to  calculate  3  for  mode  pairs  that 
have  kr  values  that  are  far  apart. 

For  adjacent  modes,  /3  can  then  be  written  in  terms  of  partial  derivatives: 


P 


dkr(n,  uj) 

__  _  _ dn _ 

did 


dkr(n,  uj) 
dn 

d2kr(n,cj) 
dn  did 


(4.10) 


4.2.2  Implicit  differentiation 

Calculating  the  partial  derivatives  in  Eq.  (4.10)  is  straightforward  if  one  has  an  ex¬ 
plicit  solution  for  kr  as  a  function  of  n  and  u.  But  explicit  solutions  for  kr  (n,  u)  can 
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only  be  obtained  for  a  few  cases  such  as  modes  in  an  ideal  waveguide  or  surface- 
trapped  modes  in  an  n2-linear  waveguide. 

Therefore,  the  present  analysis  calculates  the  partial  derivatives  in  Eq.  (4.10) 
using  implicit  differentiation.  Following  the  procedure  for  implicit  differentiation 
in  [31,  Sec.  7.2],  write  Eq.  (4.2)  as 

n -((/>  + A(/>dn  +  A(/>up  +  1)  =  0.  (4.11) 


Then  def  ne  a  function: 


f(kr,n,oj)  =  n-  (4>  +  Acpdn  +  Afup  +  1)  =  0  (4.12) 


which  allows  kr  to  be  def  ned  as  an  implicit  function  of  n  and  oj  in  regions  where 
the  derivative  df/dkr  exists  and  is  not  equal  to  zero.  The  derivative  df  /dkr  exists 
and  is  not  equal  to  zero  except  at  a  fnite  number  of  points,  such  as  kr  =  A:(0) 
and  kr  =  k(d),  as  described  in  Appendix  4. A.  So  for  our  purposes  kr  may  be 
considered  an  implicit  function  of  n  and  oj. 

The  f  rst  step  to  calculate  (3  using  implicit  differentiation  is  to  determine  the 
f  rst  order  partial  derivatives  of  kr  with  respect  to  n  and  with  respect  to  oj.  These 
can  be  obtained  using  [31,  Eq.  17], 

dkr  _  df/dn 
dn  df  /  dkr 


dkr  df  /  doj 

doj  df  /  dkr 


(4.14) 


where  df  /dx  indicates  the  result  of  differentiating  /  with  respect  to  the  explic¬ 
itly  appearing  variable  x,  holding  all  other  explicitly  appearing  variables  constant. 
More  specif  cally,  the  partial  derivatives  with  respect  to  kr  on  the  right  hand  sides 
of  Eqs.  (4.13)  and  (4.14)  should  treat  kr  as  if  it  were  not  a  function  of  n  or  oj.  See 
[31,  Secs.  7. 1-7.2]  for  details. 

To  obtain  d2kr/ (doj  dn)  for  the  numerator  ofEq.  (4.10),  the  partial  derivative 
of  (jjf  )  needs  to  be  taken  with  respect  to  oj.  This  can  be  done  as  follows: 

dk 

1.  Calculate  -  using  Eq.  (4.13). 
dn 
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d  ( dk  \ 

2.  Calculate  — —  (  — — ^  I ,  where  the  partial  derivative  with  respect  to  ui  must 
oco  \  on  J 

treat  kr  as  a  function  of  ui. 


3.  The  result  of  step  2  will  contain  terms  of  the  form  k'r(cu).  Those  terms  should 

dkr 

be  replaced  with  — —  as  calculated  by  Eq.  (4.14). 

oco 


The  results  can  the  be  inserted  into  Eq.  (4.10)  obtain  /3. 

As  will  be  discussed  in  Sec.4.6,  Eq.  (4.14)  can  be  used  to  determine  the  group 
speed  of  a  mode  at  a  particular  Ay.  Elowever  it  cannot  directly  be  used  to  obtain 
a  dispersion  curve  for  a  particular  mode,  so  the  utility  of  such  a  calculation  is  not 
immediately  apparent. 


Example 


As  an  example  of  how  to  use  the  method  described  in  this  section  to  calculate  /3, 
we  will  apply  the  method  to  an  ideal  waveguide  with  pressure-release  boundaries 
(depth  d,  sound  speed  c): 


/(Ay,  n,  to)  =  n  — 


The  application  of  Eqs.  (4.13)  and  (4.14)  yield: 


dn  d  ■  Ay 


(4.15) 


(4.16) 


dkr  co 
dco  c2kr 
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(4.17) 


Then  we  apply  the  steps  to  determine  the  mixed-partial  derivative: 


d  ( d kr\ 
duj\dn) 


TTLO  ( Ulk'r(uj )  —  kr(u)) 


Inserting  Eq.  (4.16)  and  Eq.  (4.18)  into  Eq.  (4.10)  yields: 

r2u  2  7,2 

E-  rv  f  rVj. 

n  =  =  ¥' 


(4.18) 


(4.19) 


For  modes  far  from  cut-off  kr  ~  k  leading  to  /3  «  1,  which  matches  the  result 
in  Sec.  2.3.2.  Note  that  the  estimate  of  (3  is  a  function  of  kr  instead  of  begin  a 
function  of  n  because  the  differentiation  was  performed  implicitly. 


4.3  Exploiting  dependence  on  u 

In  the  previous  section  it  was  shown  that  a  mixed  second-order  partial  derivative 
of  kr  with  respect  to  n  and  u  was  required  to  calculate  6  [Eq.  (4.10)],  implying 
that  /3  depended  on  two  independent  variables.  In  this  section  we  show  that  under 
certain  circumstances,  <j>,  A<y<m  and  A^up  depend  on  kr  and  ui  in  such  a  way  that 
the  mixed  second-order  partial  derivative  can  be  expressed  as  a  non-mixed  (regular) 
second-order  derivative,  meaning  that  fd  depends  only  on  one  independent  variable. 
It  should  be  noted  that  this  result  is  not  entirely  new:  In  [18]  Chuprov  mentions 
that  under  the  WKB  approximation,  the  functional  dependence  of  the  group  speed 
on  the  phase  speed  —  a  quantity  that  can  be  used  to  def  ne  (3  —  does  not  depend 
on  the  frequency  under  certain  circumstances. 

We  begin  by  noting  that  a  plot  of  4>(kr,uj)  will  look  like  a  scaled  down  version 
of  a  plot  of  <j)(kr,  oj  +  Ac u);  both  axes  shrink  uniformly.  Furthermore,  the  shape  of 
A^dn  =  (l>b  also  scales  with  frequency,  but  the  scaling  only  occurs  on  the  to  axis. 
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This  scaling  suggests  that  the  shape  of  <f>(kr,Lj)  and  (Ay,  oj)  at  one  frequency 
as  a  function  of  kr  provides  suff  cient  information  to  take  partial  derivatives  with 
respect  to  either  kr  or  oj  at  any  frequency.  We  now  show  that  0(kr.oj)  can  be 
written  as  oj  ■  0(1,  kr/oj)  and  A0}n(Ay,  oj)  can  be  written  as  A<0 n(kr/oj). 


First,  factor  an  oj  out  of  kz(z): 


(4.20) 


(4.21) 


Note  that  (^)  is  the  inverse  of  the  phase  speed,  which  is  known  as  the  phase 
slowness  Sp. 

Sp  =  -  (4.22) 

OJ 

The  remainder  of  this  chapter  will  denote  (^)  by  Sp  to  simplify  notation,  but  it  is 
important  to  keep  in  mind  the  Sp  is  a  function  of  u  when  taking  derivatives  with 
respect  to  oj. 

Next  we  show  that  0^,  the  phase  of  the  ref  ection  coeff  cient  between  two- 
homogeneous  half-spaces,  and  thus  <j>^  can  be  written  as  a  function  of  Sp.  Fol¬ 
lowing  [33,  Eq.  1.45  and  Eq.  2.125], 

Z>  —  /f 

Ref  ection  coeff  cient  =  R  = -  (4.23) 

Zi  +  Z2 
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where 


Pi^i 
sin  Oi 
kpiCi 
kz 

(' u/Ci)pjCj 
yJ{u/Ci)2  -  k2 

_ Pi _ 

\/{l/ci)2  ~  (kr/u)2 
Pi 

sj{l/ci)2  -  S 2 


(4.24) 


Then  R  can  be  written  as  a  function  of  Sp,  and  so  can  Aqyn. 

To  see  that  A0up  (and  A^dn  in  the  case  of  a  vacuum  bottom)  can  be  written  as 
function  of  Sp,  note  that  the  value  of  A(pup(kr,  ui)  as  defned  in  Sec.  4.2  depends 
only  on  the  point  at  which  kr  =  or  equivalently  at  Sp  = 

Having  written  0,  A0up,  and  A(f>&a  as  function  of  ( kr/u ;),  the  waveguide  in¬ 
variant  j3  can  be  calculated  by  writing  Eq.  (4.12)  as 


0  =  f(kr,n,u) 

=  n-  ■  (f>\  +  A^dn  +  A<^up  ^  (4-25) 

and  then  following  the  steps  in  Sec.  4.2.2  to  calculate  /?  from  f(kr,n,u).  The 
resulting  expressions  are  lengthy  and  diff  cult  to  manipulate  by  hand,  so  the  author 
used  the  computer  algebra  system  Mathematica.  See  Appendix  4.B  for  details.  The 
f  nal  result  can  be  simplif  ed  down  to: 


g(S  u)  = _ Ml  (Sp)  +  (^))2 _ 

A^dn  ( SP  f  +  wA (t)'An  (Sp)  </>[  (Sp)  +  w</>  1  (Sp)  (Ad>"n  (Sp)  +  U(t> '/  (Sp)) 

(4.26) 

where  the  '  denotes  a  derivative  with  respect  to  the  argument  Sp  =  Note  that 
A(j)'ln  is  related  to  the  “ray  displacement”  discussed  in  [1 1,  Sec.  4.4], 

When  Sp  >  1  /c(d)  (i.e,  when  the  mode  has  a  lower  turning  point)  or  if  Ac^n 
represents  a  vacuum  bottom,  then  the  derivative  of  A^  with  respect  to  Sp  is  zero 


everywhere  except  possibly  at  one  point,  and  Eq.  (4.26)  can  be  reduced  to 


P(Sp)  = 


(SP)2 

d>i(sp)d/i(spy 


(4.27) 


This  expression  is  remarkably  simple  given  the  complexity  of  its  derivation,  which 
is  shown  in  Appendix  4.B.  When  it’s  valid,  Eq.  (4.27)  provides  an  direct  rela¬ 
tionship  between  the  SSP  and  the  value  of  /3  at  all  horizontal  wavenumbers  and 
all  frequencies  though  a  single  function  of  a  single  variable:  (j>i(Sp).  It  should  be 
noted  that  Eq.  4.27  is  equivalent  to  [14,  Eq.  24]. 


4.4  Interpretation  and  implementation 


Eq.  (4.26)  reveals  that  the  shape  of  <f>i(Sp)  and  <f>4n(Sp)  provides  suffcient  infor¬ 
mation  to  calculate  the  value  of  j3  for  any  kr  and  u.  Eq.  (4.27)  reveals  that  when 
Sp  >  1/ c(d )  or  when  the  bottom  is  approximated  as  being  a  vacuum,  B  depends 
only  the  ratio  of  kr  to  u>,  as  opposed  to  depending  on  each  variable  independently. 

To  relate  the  value  of  f3  to  the  SSP  c(z)  and  Sp,  some  properties  of  B\ (Sp)  will 
now  be  discussed.  Let  cmm  and  cmax  be  the  minimum  and  maximum  sound  speed 
in  the  water  column,  not  including  the  seaf  oor.  cf>i  (Sp)  can  be  broken  up  into  two 
regions,  the  frst  being  (0  <  Sp  <  l/cmax)  which  corresponds  to  SRBR  modes  (see 
Sec.  2.2.2).  And  the  second  being  (l/cmax  <  Sp  <  1/cmin)  which  corresponds  to 
non-SRBR  modes.  In  both  regions  B]  (Sp)  is  a  positive  and  monotonically  decreas¬ 
ing  function. 

Modes  with  horizontal  wavenumbers  in  the  f  rst  region  do  not  have  any  turning 
points  and  “feel”  the  top  and  bottom  boundary  of  the  waveguide.  In  this  region 
B\  (Sp)  and  its  derivatives  can  be  approximated  to  arbitrary  accuracy  by  breaking 
the  water-column  up  into  several  small  iso-velocity  segments,  or  equivalently,  ap- 
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proximating  Eq.  (4.21)  with  Riemann  Sum: 


MSp) 


MSP) 

ti(Sp) 


(4.28a) 

(4.28b) 


(4.28c) 


where  <j>\ ,i(Sp)  represents  the  contribution  from  the  i-th  segment,  and  di  is  the  depth 
of  the  i-th  segment.  For  ( Sp  <  l/cmax),  every  term  in  Eq.  (4.28a)  is  a  positive  but 
monotonically  decreasing  function  with  negative  concavity,  so  4>"(SP)  <  0,  as  can 
be  verifed  by  inspection  of  Eq.  (4.28c).  Inserting  the  sign  of  4>\(SP)  and  4>"(SP) 
into  Eq.  (4.27)  reveals  that  (3  >  0  for  all  SRBR  modes  in  a  waveguide  with  a 
pressure-release  boundaries,  regardless  of  the  SSR 

For  non-SRBR  modes  (l/cmax  <  Sp  <  l/cmjn),  Eqs.  (4.28b)  and  (4.28c)  are 
incorrect  because  they  do  not  take  into  account  that  the  limits  of  the  def  nite  integral 
of  kz(z)  in  Eq.  (4.21)  —  the  turning  points  —  are  themselves  functions  of  Sp.  3 
Therefore  in  SSP  segments  where  a  turning  point  exists,  c(z)  must  be  represented 
with  a  segment  that  is  not  iso-velocity  so  that  the  derivative  with  respect  to  Sp  can 
take  into  account  that  the  limits  of  the  integral  in  Eq.  (4.21)  are  functions  of  Sp.  For 
this  we  use  n2- linear  segments  ([33,  Sec.  2.5.1]  and  [5]),  which  can  be  analytically 
integrated  and  are  very  close  to  linear  for  the  SSPs  typically  encountered  in  ocean 
acoustics.  A  segment  with  an  n2- linear  SSP  is  one  of  the  form 

°i(z)  =  7  1  ,  (4-29) 

VciiZ  +  bi 


where  a*  and  bi  are  chosen  to  interpolate  the  SSP.  Each  segment  can  have  its  own 
coordinate  system,  so  without  loss  of  generality  let  <  0. 


3See  “Leibniz  Integral  Rule”  in  any  introductory  Calculus  textbook.  Note  that  although  the 
value  of  <f>i(Sp)  is  not  affected  by  the  fact  that  the  integration  limits  are  functions  of  Sp  because 
kz(z,  Sp)  —  0  when  ( z ,  Sp)  corresponds  to  a  turning  point,  the  value  of  4>'{(SP)  is  affected. 
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Inserting  Eq.  (4.29)  into  Eq.  (4.21)  yields 


MSpi  =  ‘ ~ ^) dz 

=  db  + aidi  ~  s^3/2  ~  ( bi  ~  5p)3/2)  •  (4-30) 

For  a  segment  that  contains  a  turning  point,  ^ Jbi  +  aidi  —  5|  =  0  at  that  turning 
point  and  thus 

MsP)  =  ^{bi-slf/2 .  (4.31) 

To  represent  an  arbitrary  SSP,  one  can  use  any  combination  of  segments  that  are 
iso-velocity  or  n2- linear  (or  any  c(z)  for  which  there  is  a  closed  form  expression  of 
Eq.  (4.21)),  as  long  as  the  segments  that  contain  turning  points  are  not  iso-velocity. 

4>i (Sp),  <j>i( Sp ),  4>1(Sp),  which  are  needed  to  calculate  /?,  can  then  be  calcu¬ 
lated  quasi-analytically  in  the  sense  that  the  derivatives  are  taken  analytically  for 
each  term  of  the  summation  that  represents  c(z),  but  the  summations  themselves 
may  be  computed  numerically  if  there  are  too  many  terms  to  handle  analytically. 
In  the  case  of  a  waveguide  with  a  bottom  half-space,  the  function  4>b{Sp )  is  also  re¬ 
quired  to  calculate  3  and  its  derivatives  can  calculated  analytically  from  Eq.  (4.23). 

Many  waveguides  can  be  suff  ciently  well  represented  with  just  a  few  seg¬ 
ments  and  a  bottom  half-space,  and  thus  j3  can  be  computed  purely  analytically 
in  those  cases.  In  all  waveguides,  regardless  of  how  complicated  the  SSP  is,  4>i  has 
a  straightforward  meaning  in  the  sense  that  it  is  easy  to  qualitatively  understand 
how  the  details  of  the  SSP  affect  4>\ . 


4.5  Examples 

This  section  applies  the  methods  described  in  Secs.  4.2  -  4.4  to  calculating  j3  in 
waveguides  with  several  commonly  analyzed  SSPs:  iso-velocity,  a  single  n2- linear 
segment,  and  a  deep-water  Munk  SSP.  For  each  SSP,  3  is  plotted  for  a  vacuum 
bottom  and  for  a  fuid  halfspace  bottom.  In  addition  k(z),  3\  and  dyn  are  plotted 
on  an  x-axis  common  with  the  plots  of  3  so  that  their  relationship  to  3  can  be 
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easily  seen.  In  all  cases  that  have  a  fuid  bottom  halfspace,  Cbottom  =  1650m/s, 

Pbottom  —  1.5pwater- 

For  comparison,  the  normal  mode  program  Kraken  [48]  was  also  used  to  cal¬ 
culate  the  value  of  3  using  the  same  environmental  parameters.  3  was  calculated 
from  Kraken  using  Eq.  2.27,  using  only  adjacent  mode  numbers. 

All  calculations  in  Kraken  used  an  acoustic  source  frequency  of  250  Hz  except 
for  the  deep-water  Munk  case  which  used  a  frequency  of  75  Hz.  As  the  frequency 
increases,  the  agreement  between  Kraken  and  the  methods  described  in  this  chapter 
tends  to  improve,  as  would  be  expected  because  the  accuracy  of  the  WKB  approx¬ 
imation  increases  with  frequency  as  does  the  accuracy  of  approximating  the  f  nite 
difference  in  mode  order  with  a  continuous  derivative.  This  agreement  provides  a 
sanity  check  on  the  validity  of  the  present  analysis. 

4.5.1  Iso-velocity  SSP:  ideal  and  Pekeris  waveguides 

An  iso-velocity  SSP  with  a  vacuum  bottom  is  known  as  an  ideal  waveguide.  Al¬ 
though  3  for  an  ideal  waveguide  can  be  calculated  analytically  using  Eq.  (2.15) 
because  explicit  expressions  for  kr(m,u )  can  be  obtained  (See  Sec.  2.3.2),  the 
result  is  reproduced  here  for  illustrative  purposes.  For  an  ideal  waveguide 

Msp)  =  ^\J^-S^  (4.32) 

Calculating  (£>[( SP )  and  4>"(SP)  and  inserting  into  Eq.  (4.27)  yields: 

3  =  c2Sp.  (4.33) 

As  Sp  -»  l/c,  which  happens  as  modes  move  away  from  being  cut  off,  3  — >•  1,  as 
shown  in  Fig.  4.4(c). 

An  iso-velocity  SSP  with  a  bottom  f  uid  halfspace  is  known  as  a  Pekeris  waveg¬ 
uide.  4>i(Sp)  is  the  same  as  that  for  an  ideal  waveguide  but  the  derivatives  of 
4>dn{Sp)  are  non-zero,  so  Eq.  (4.26)  must  be  used  instead  of  Eq.  (4.27).  4>dn(Sp) 
for  a  bottom  that  is  denser  and  has  a  faster  c  than  the  water  just  above  it  always 
has  the  same  qualitative  shape,  shown  in  Fig.  4.4(b).  <3An(Sp)  and  Kn (sp)  can  be 
calculated  analytically  in  a  straightforward  manner,  although  the  result  is  a  lengthy 
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expression.  Eq.  (4.26)  can  then  be  used  to  obtain  an  analytic  estimate  of  (3  as  a 
function  of  Sp.  In  the  case  of  a  Pekeris  waveguide,  the  WKB  approximation  is  the 
exact  solution,  so  the  only  approximation  made  in  calculating  (3  is  approximating 
the  f  nite  difference  in  mode  number  with  a  continuous  derivative.  The  reason  one 
can  obtain  an  analytic  estimate  of  /3  in  a  Pekeris  waveguide  despite  not  being  able 
to  analytically  solve  for  kr(m,  co)  is  that  the  resulting  /3  is  a  continuous  function 
of  Sp  instead  of  being  an  explicit  function  of  mode  number.  Figure  4.4(c)  shows  a 
plot  of  f3  for  a  Pekeris  waveguide. 
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WKB-based  approximation  (vacuum  bottom) 
Kraken  (vacuum  bottom) 

WKB-based  approximation  (fluid  bottom) 
Kraken  (fluid  bottom) 


Figure  4.4:  The  waveguide  invariant  j3  for  an  ideal  and  a  Pekeris  waveguide,  and 
the  functions  used  to  calculate  (3,  all  plotted  on  a  common  x-axis  scale  (the  phase 
slowness  Sp  normalized  by  the  minimum  sound  speed  in  the  water  column),  (a) 
k(z)  =  which  in  this  case  is  a  constant  and  equal  to  1  because  of  the  nor¬ 
malized  scale,  (b)  The  (j>i(Sp)  associated  with  k(z)  (black),  and  0dn (Sp)  for  the 
bottom  fuid  halfspace  (gray),  (c)  (3  versus  phase  slowness.  For  the  ideal  waveg¬ 
uide  (vacuum  bottom)  the  WKB-based  approximation  of  (3  closely  matches  the 
value  calculated  from  Kraken.  For  the  Pekeris  waveguide  (fuid  bottom)  there  is 
also  good  agreement  between  Kraken  and  the  WKB-based  approximation.  For  the 
f  uid  bottom,  the  WKB-based  approximation  of  f3  is  unbounded  at  the  value  of  Sp 
corresponding  to  the  critical  angle  because  the  derivative  of  4>&n (Sp)  is  infnite. 
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4.5.2  n2  -linear  waveguide 

The  n2 -linear  waveguide  used  in  this  section  has  c(0)  =  1500  m/s  and  c(d )  =  1550 
m/s,  and  d  =  100  m.  The  k(z )  associated  with  c(z )  is  shown  in  Fig.  4.5(a). 


Surface-trapped  modes 


Surface-trapped  modes  are  modes  that  have  Sp  >  1  /c(d)  in  an  upward  refract¬ 
ing  SSP.  Using  Eqs.  (4.31)  and  (4.25),  one  can  obtain  an  explicit  solution  for 
kr(m cu),  use  a  Taylor  series  expansion  similar  to  that  in  [1 1,  Sec.  6.7.2]  and  then 
use  Eq.  (2.15)  to  see  that  ft  ~  —3  (or  modes  far  from  cut-off.  Alternatively  one  can 
use  the  method  in  Sec.  4.3.  Using  Eq.  (4.31)  to  represent  and  inserting  (ft  into 


Eq.  (4.27)  yields: 


P 


3Sp 

b-2Sj' 


(4.34) 


Surface-trapped  modes  far  from  cutoff  have  a  phase  slowness  of  Sp  «  1  / c(0)  = 
Vb.  Inserting  that  into  Eq.  (4.34)  results  in  ft  =  —  3. 

For  Sp  >  1  /c(d)  the  concavity  of  (ft  is  positive,  as  shown  in  Fig.  4.5(b).  This 
is  because  as  Sp  decreases,  the  integrand  increases  in  value  and  the  limits  of  the 
integral  move  further  apart.  The  result  is  that  ft  is  a  negative  number  for  Sp  > 
1  /c(d)  as  can  be  seen  in  Fig.  4.5(c). 

This  analysis  for  surface-trapped  modes  also  applies  to  waterborne  modes  in  a 
two-segment  ducted  waveguide  where  each  side  of  the  duct  has  an  n2- linear  SSP, 
because  such  modes  have  a  (ft  equivalent  to  that  of  the  surface-trapped  modes 
discussed  above,  and  so  those  modes  will  also  have  ft  «  —  3  [26]. 


SRBR  modes 

SRBR  modes  in  this  upward  refracting  environment  have  Sp  >  l/c(d).  Based  on 
the  qualitative  argument  in  Fig.  4.1,  one  would  expect  that  ft  for  the  SRBR  modes 
should  be  similar  to  ft  in  an  iso-velocity  waveguide  when  (k( 0)  —  k(d))/(k( 0)  — 
kr  )  <C  1.  An  explicit  analytic  solution  for  kr(m,u )  of  these  modes  is  either  in¬ 
tractable  or  does  not  exist.  Flowever  ft  can  be  calculated  with  fairly  simple  expres¬ 
sions  using  Eqs.  (4.30)  and  (4.27).  ft  for  a  vacuum  and  a  f  uid  halfspace  bottom  are 
shown  in  Fig.  4.5(c). 
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-  WKB-based  approximation  (vacuum  bottom) 

-  -  -  Kraken  (vacuum  bottom) 

WKB-based  approximation  (fluid  bottom) 

Kraken  (fluid  bottom) 

Figure  4.5:  The  waveguide  invariant  B  for  an  n2- linear  waveguide,  and  the  func¬ 
tions  used  to  calculate  B  plotted  on  a  common  x-axis  (the  phase  slowness  Sv  nor¬ 
malized  by  the  minimum  sound  speed  in  the  water  column),  (a)  Normalized  k(z). 
(b)  The  4>i  (Sp)  associated  with  k(z)  (black),  and  4>dn{Sp )  for  the  bottom  fuid  half¬ 
space  (gray).  Note  that  the  concavity  of  4>i(Sp)  changes  when  Sp  =  l/c(d),  which 
results  in  a  sign  change  of  (3  at  that  point,  (c)  B  for  the  vacuum  and  f  uid  bottom 
n2- linear  waveguides.  The  WKB-based  approximation  is  in  good  agreement  with 
Kraken  everywhere  except  in  the  region  near  Sp  =  1  /c(d)  =  0.974  •  cm;n  which  is 
where  the  mode’s  turning  point  is  close  to  the  bottom  and  the  WKB  approximation 
is  inaccurate.  For  Sp  >  l/c(d)  the  modes  are  surface-trapped,  B\  has  positive  con¬ 
cavity,  and  B  ~  —3.  For  Sp  <  l/c(d),  which  corresponds  to  bottom-interacting 
modes,  (3  «  1.  The  effect  of  including  a  bottom  f  uid  halfspace  instead  of  a  vacuum 
is  small,  except  for  modes  close  to  cutoff  (modes  with  kr  values  near  the  critical 
angle  kr  value)  because  <fi'dn(Sp)  becomes  unbounded  at  that  point. 
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4.5.3  Quasi-analytic  results  for  a  deep-water  Munk  profile 

A  deep-water  waveguide  with  a  Munk  SSP  [33,  Sec.  5.6]  can  be  approximated  with 
n2 -linear  segments.  Fig.  4.6(c)  shows  3  for  5000  m  deep  waveguide  with  a  Munk 
SSP,  for  both  a  vacuum  and  a  f  uid  halfspace  bottom.  The  (l>\  used  for  calculating 
f3  was  computed  with  Eq.  (4.30)  using  1000  segments  to  represent  c(z). 
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- WKB-based  approximation  (vacuum  bottom) 

- Kraken  (vacuum  bottom) 

WKB-based  approximation  (fluid  bottom) 
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Figure  4.6:  The  waveguide  invariant  (3  for  a  deep-water  Munk  SSP,  and  the  func¬ 
tions  used  to  calculate  (3  plotted  on  a  common  x-axis  (the  phase  slowness  Sp 
normalized  by  the  minimum  sound  speed  in  the  water  column),  (a)  Normalized 
k(z).  (b)  The  <j>i(Sp)  associated  with  k(z')  (black),  and  4>dn{Sp)  for  the  bottom 
f  uid  halfspace  (gray),  (c)  (3  for  the  Munk  SSP  with  vacuum  and  f  uid  bottoms. 
In  general  there  is  good  agreement  between  the  WKB-based  approximation  and 
Kraken  except  at  Sp  values  where  the  mode’s  turning  point  is  near  a  boundary  and 
the  WKB  approximation  is  inaccurate.  In  regions  where  modes  are  non-SRBR 
(Sp  ■  Cmin  >  0.968),  the  concavity  of  4>i  (Sp)  is  negative  and  so  (3  is  positive.  For 
Sp  values  corresponding  to  top-  and  bottom-interacting  modes,  3  is  between  about 
0.5  and  1.  The  effect  of  including  bottom  fuid  halfspace  instead  of  a  vacuum  is 
negligible,  which  is  why  no  gray  lines  can  be  seen. 
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4.5.4  Realistic  shallow  water  waveguide 


Section  3.4.2  presented  some  experimentally  measured  waveguide  invariant  stria- 
tions  from  a  shallow  water  waveguide  (see  that  section  for  details).  In  this  subsec¬ 
tion,  we  calculate  the  value  of  (3  for  that  environment  (i.e.  for  the  SSP  in  Fig.  3.5, 
assuming  Cb0ttom  =  1500  m/s).  It  can  be  seen  that  3  has  a  value  slightly  less  than 
one  for  the  SRBR  modes,  and  that  (3  varies  greatly  for  non-SRBR  modes. 
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- WKB-based  approximation  (vacuum  bottom) 

- Kraken  (vacuum  bottom) 

- WKB-based  approximation  (fluid  bottom) 

Kraken  (fluid  bottom) 


Figure  4.7:  The  waveguide  invariant  /3  for  a  realistic  shallow  water  waveguide,  and 
the  functions  used  to  calculate  f3  plotted  on  a  common  x-axis  (the  phase  slowness 
Sp  normalized  by  the  minimum  sound  speed  in  the  water  column).  The  SSP  for 
this  fgure  comes  from  Fig.  3.5  (a)  Normalized  k(z).  (b)  The  0\  (Sp)  associated 
with  k(z)  (black),  and  (t>dn(SP)  for  the  bottom  fuid  halfspace  (gray),  (c)  3  for  the 
realistic  shallow  water  SSP  with  vacuum  and  f  uid  bottoms.  In  regions  where  modes 
are  SRBR  ( Sp  ■  rmm  <  0.98),  3  is  slightly  less  than  1.  For  Sp  values  corresponding 
SRBR  modes  (Sp  ■  cmm  >  0.98),  the  value  of  3  varies  greatly. 
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4.6  Calculating  (3  for  non-adjacent  modes 


The  def  nition  of  /3  in  its  most  general  form,  Eq.  (2. 15),  is  for  a  pair  of  modes  whose 
kr  values  are  not  necessarily  adjacent  to  one  another.  The  methods  for  calculating 
3  in  Secs.  4.2  -  4.5  are  valid  only  for  mode  pairs  that  have  kr  values  that  are  close 
enough  together  so  that  the  curve  (l>\  (Sp  =  kr/u)  is  approximately  linear  between 
the  two  modes’  Sp  values.  This  restriction  on  the  validity  is  because  the  fnite 
difference  of  the  modes’  kr  values  was  approximated  to  be  inf  nitesimally  small 
(i.e.,  a  derivative).  Because  of  this,  3  depended  only  on  a  single  kr  value.  The 
justif  cation  for  such  an  approach  was  discussed  in  the  beginning  of  Sec.  4.2.1. 

However,  some  components  of  the  acoustic  intensity  f  eld  I (?’,  u)  will  be  the 
result  of  mode  pairs  that  have  Sp  values  that  are  far  apart  (non-adjacent).  When  cal¬ 
culating  3  for  such  mode  pairs,  the  f  nite  difference  between  the  modes’  kr  values 
cannot  be  approximated  by  a  derivative  as  was  done  in  Secs.  4.2  -  4.5.  However 
some  of  the  analysis  from  Secs.  4.2  and  4.3  is  still  valid  for  a  pair  of  non-adjacent 
modes,  and  can  be  applied  to  calculate  3  for  non-adjacent  modes.  The  resulting 
3  will  depend  explicitly  on  the  two  modes’  Sp  values,  and  implicitly  on  the  two 
modes’  modal  numbers. 

Using  Eq.  (2.15)  as  a  def  nition  for  3,  Eq.  (4.25)  can  be  inserted  into  Eq.  (4.14) 
to  implicitly  the  calculate  derivatives  with  respect  to  u  in  the  denominator  Eq.  (2.1 5). 


9K  =  <I>1  (fe)  ~  M  (£)  ~  (M  (4  ^ 

If  the  bottom  f  uid  halfspace  is  approximated  as  being  a  vacuum,  then  <j)'dn  =  0  and 


dkr  _  (f>i  {Sp)  +  Sp(j> 1  (Sp) 

du  <?!>i  (sp) 

-  ( Q 

V  p  #  (Sp) 


(4.36) 


Note  that  this  is  the  group  slowness.  It  is  worth  mentioning  that  Eq.  (4.36)  is 
identical  to  [45,  Eq.  2.14.12],  where  it  appears  in  the  context  of  ray  theory  and  Sp 
was  def  ned  as  the  slowness  of  a  ray  at  its  turning  point,  as  opposed  to  the  def  nition 
used  here  which  is  in  terms  of  the  modal  phase  slowness  and  thus  explicitly  takes 
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into  account  frequency  dependence. 

Eq.  (4.35)  or  (4.36)  can  then  be  inserted  into  Eq.  (2.15)  to  calculate  p  for  a  pair 
of  modes  with  Sp  values  that  are  not  adjacent.  When  the  bottom  is  a  vacuum: 


P (.kr, mi  kr,h  ^ ) 


1  ( krm  fc r,l ) 

UJ  d(kr,m  ^r,i) 
di o 

—  ( kr{) 

(C  _  C  _ 

[bp’m  bp'l)+  ^i(sP,0 

_ {Sp,m  ~  Spj) _ 

/  Q  _  Q  \  |  f  <t>l  {Sp,l)  _  01 (Sp  ,m )  \ 
-  3p,J)  -E  ^0^5  -  0i(Sp,m)  J 


(4.37) 


P  (Sp.m,  Sp,i)  can  then  be  plotted  as  a  function  of  the  two  modes’  phase  slow¬ 
nesses.  Although  the  specif  c  SPjTn  and  Spj  values  that  correspond  to  propagating 
modes  are  unknown  unless  one  explicitly  solves  for  them,  one  can  determine  the 
regions  of  P  (<SP)rn,  Spj)  that  correspond  to  SRBR  modes  and  non-SRBR  modes. 

Fig.  4.8  shows  P  {SPtm,  Spj)  for  the  same  n2-linear  waveguide  with  a  vacuum 
bottom  used  in  the  previous  section,  calculated  using  Eq.  (4.37)  and  using  Kraken 
for  comparison.  A  fuid  bottom  could  be  taken  into  account  inserting  Eq.  (4.35) 
into  Eq.  (2.15). 


4.7  Chapter  summary 

Knowing  the  value  of  13  is  critical  for  many  applications  of  the  waveguide  invari¬ 
ant,  such  as  the  passive  range  estimation  presented  in  Ch.  3.  Under  the  WKB 
approximation,  (3  can  be  determined  solely  from  the  phase  of  the  bottom  Rayleigh 
ref  ection  coeff  cient  as  a  function  of  the  phase  slowness  and  from  the  integrated 
phase  of  kz  in  the  water  column  as  a  function  of  the  phase  slowness,  (!>\  (Sp).  If 
the  waveguide  has  a  bottom  f  uid  halfspace,  P  depends  on  kr  and  ui  independently 
and  can  be  calculated  using  Eq.  (4.26).  The  use  of  implicit  differentiation  allows 
for  analytic  estimates  of  P  in  waveguides  whose  horizontal  wavenumbers  cannot 
be  explicitly  solved  for  analytically.  When  the  bottom  is  approximated  as  being 
a  vacuum,  f3  depends  only  on  the  ratio  of  kr  to  u  (phase  slowness),  has  a  simple 
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Figure  4.8:  Plots  of  the  waveguide  invariant  j3  for  non-adjacent  modes  versus  nor¬ 
malized  phase  slowness  (Sp  ■  cmm)  for  each  mode  (S[hrn  and  Spj).  These  plots  are 
from  the  same  n2 -linear  SSP  used  for  Fig.  4.5,  with  a  vacuum  bottom  and  were 
calculated  from  (a)  Kraken  and  (b)  Eq.  (4.37).  The  white  areas  along  the  diagonal 
correspond  to  m  =  (in  Eq.  (2.15)  and  are  thus  undefned.  A  diagonal  “slice”  just 
off  the  line  Sp.m  =  Spj  corresponds  to  adjacent  modes,  and  is  the  same  as  what  is 
plotted  in  Fig.  4.5(c)  for  the  vacuum  bottom.  For  this  SSP  f3  «  1  for  all  mode  pairs 
where  both  modes  are  SRBR,  and  3  ~  —3  lor  all  mode  pairs  where  both  modes 
are  non-SRBR.  In  regions  where  one  mode  is  SRBR  and  the  other  is  non-SRBR,  3 
usually  greater  than  1 . 
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relationship  to  the  integrated  vertical  phase  <f>i(Sp),  and  can  be  calculated  using 
Eq.  (4.27)  for  adjacent  modes  and  Eq.  (4.37)  for  non-adjacent  modes. 

The  relationship  between  c(z),  <f>i(Sp),  and  /3  provides  a  means  for  understand¬ 
ing  how  the  SSP  affects  the  value  of  f3.  A  qualitative  argument  was  given  which 
suggests  that  (3  ~  1  for  SRBR  modes  in  waveguides  where  the  seaf  oor  is  modeled 
as  a  vacuum,  regardless  of  the  SSR  It  was  also  shown  that  f3  for  non-SRBR  modes 
can  take  on  a  wide  range  of  values,  which  depend  strongly  on  the  details  of  the  SSP. 

The  present  analysis  discussed  how  to  calculate  (3  for  a  single  pair  of  modes. 
But,  as  can  be  seen  in  Eq.  (2.9),  the  acoustic  intensity  is  a  summation  of  terms 
from  all  mode  pairs.  Each  term  of  the  summation  will  contribute  to  the  total  inten¬ 
sity  (although  a  few  groups  of  modes  sometimes  dominate  the  summation).  The 
relative  amplitude  of  each  term  of  the  summation  depends  on  the  relative  modal 
amplitudes,  which  in  turn  depend  primarily  on  the  source  and  receiver  depths  rel¬ 
ative  to  features  in  the  sound  speed  prof  le.  This  depth-dependence  is  discussed 
brief  y  in  [51]  and  [52].  For  example,  if  both  the  source  and  receiver  are  near  a  lo¬ 
cal  sound  speed  minimum,  the  local  modal  amplitudes  of  the  modes  trapped  in  that 
sound  speed  minimum  will  likely  be  large,  and  one  will  observe  some  striations 
according  to  the  value  of  f3  for  mode  pairs  with  Sp  «  1/cmin.  The  present  analysis 
provides  a  rigorous  framework  for  analyzing  how  the  observed  striation  slopes  in 
an  acoustic  intensity  plot  depend  on  the  source  and  receiver  depths  and  the  SSR 

As  with  many  ocean  acoustic  phenomena,  the  waveguide  invariant  can  be  un¬ 
derstood  from  the  perspective  of  normal  modes  or  ray  theory.  The  insights  ob¬ 
tained  by  the  present  normal  mode  analysis  complement  the  ones  obtained  by  the 
ray-theoretic  analysis  presented  in  [24,  22,  14].  In  fact,  because  ray  theory  can  be 
derived  by  taking  the  WKB  approximation  to  the  limit  of  inf  nite  frequency,  the 
results  from  the  present  analysis  as  ui  — >•  oo  should  match  the  results  derived  in 
[24,  22].  Reference  [14]  discusses  the  relationship  between  normal  modes  and  ray 
theory  in  the  context  of  the  waveguide  invariant. 

Because  the  striations  from  the  waveguide  invariant  can  be  seen  in  the  beam- 
formed  output  from  a  horizontal  array,  as  will  be  discussed  in  Ch.  5,  one  possible 
application  of  the  present  analysis  would  be  to  use  a  horizontal  array  to  observe 
striations  corresponding  to  specif  c  range  of  kr  values,  and  thus  a  specif  c  range  of 
/3  values.  This  idea  is  explored  in  Ch.  5.  Knowledge  of  which  kr  values  correspond 
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to  which  f3  values  may  facilitate  passive  range  estimation  and  other  applications  of 
the  waveguide  invariant,  especially  in  low  SNR  environments. 

4.A  Chapter  Appendix:  Justification  for  treating  Ay  as  an 
implicit  function 

Ay  can  be  treated  as  an  implicit  function  of  n  and  u  in  regions  where  the  derivative 
df/dkr  exists  and  is  not  equal  to  zero  [31,  Sec.  7.2].  The  analysis  in  this  chapter 
is  only  interested  in  the  region 

OJ  i o 

<  Ay  <  Amax  =  (4.38) 

^bottom  ^min 

where  Cbottom  is  the  speed  of  sound  in  the  bottom  f  uid  halfspace  (use  oo  for  a 
vacuum),  and  cmm  is  the  minimum  sound  speed  in  the  water  column.  To  determine 
when  df  /dkr  exists  and  is  not  equal  to  zero,  use  Eq.  (4.12)  to  write: 

df  _  ( d <j>(kr,u)  dAcf)dn  dAfup\ 

dkr  \  dkr  dkr  dkr  ) 

Inspection  ofEq.  (4.4)  reveals  that  df(kr,uj)  /dkr  exists  everywhere  inside  the 
region  of  interest  except  perhaps  at  Ay  values  corresponding  to  discontinuities  in 
k(z )  =  u>/c(z)  and  at  Ay  =  A( 0)  and  Ay  =  k(d).  dAcj)Up/dkr  equals  zero  every¬ 
where  except  at  Ay  =  A( 0),  where  the  derivative  does  not  exist.  dA<fdn/dkr  can  be 
calculated  analytically  and  is  smooth,  continuous,  and  monotonically  decreasing 
inside  the  region  of  interest. 

Therefore,  Ay  can  be  def  ned  as  an  implicit  function  of  n  and  uj  in  the  region 
of  interest  excluding  the  points  Ay  =  /;•(()),  Ay  =  kid).  Ay  values  corresponding  to 
discontinuities  in  k(z),  and  possibly  one  more  location  where  df  /dkr  =  0. 

All  of  the  aforementioned  discontinuities  are  artifacts  of  the  WKB  approxi¬ 
mation  used  in  this  chapter.  In  a  full-wave  solution  of  the  depth-separated  wave 
equation,  the  effect  of  the  surface  or  sea  f  oor  on  the  horizontal  wavenumber  would 
smoothly  increase  as  the  turning  point  of  the  mode  approaches  the  waveguide 
boundary.  Future  research  could  attempt  to  extend  the  current  analysis  to  a  WKB 
approximation  that  includes  the  evanescent  Airy  tails  of  the  non-SRBR  modes. 
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4.B  Chapter  Appendix:  Derivation  of  expression  for  (3 


Equation  (4.26),  which  includes  the  effects  of  a  bottom  f  uid  halfspace,  was  derived 
with  the  assistance  of  the  computer  algebra  system  Mathematica.  As  described  in 
Sec.  4.3,  when  the  bottom  halfspace  is  a  vacuum,  Eq.  (4.26)  reduces  to  Eq.  (4.27). 
Alternatively,  Eq.  (4.27)  can  be  derived  by  hand  if  one  begins  the  derivation  assum¬ 
ing  that  the  bottom  halfspace  is  a  vacuum.  That  derivation  is  performed  here. 


Replace  the  k'r  (ut)  term  in  Eq  .(4.44)  with  Eq.  (4.43)  and  simplify  the  result: 


d  ( dkr\ 

du>  \  dn  J 


(4.45) 


Insert  Eqs.  (4.42)  and  (4.45)  into  Eq.  (4.10),  and  write  ^  as  Sp: 


P(SP)  = 


#  (^p)2 

<h  (sp)<ft(sp) 


(4.46) 


which  is  the  same  as  Eq.  (4.27). 

For  a  waveguide  with  a  fuid  bottom  halfspace,  the  expression  for  n  is  that 
given  in  Eq.  (4.25)  instead  of  the  one  given  in  Eq.  (4.40).  The  same  mathematical 
operations  are  then  be  applied,  but  the  intermediate  expressions  are  lengthy  and 
diff  cult  to  manipulate  by  hand.  The  following  lines  of  Mathematica  code  were 
used  to  produce  the  result  in  Eq.  (4.26). 

ff  =  n  -  (w  phil [kr/w]  +phidn [kr/w]  +  1) 
dkdn  =  -D [f f , n] /D [ff , kr]  (*  dkr  /  dn  *) 
dkdw  =  -D [ f f , w] /D [ f f , kr ]  (*  dkr  /  dw  *) 
stepl  =  dkdn/ . kr->kr [ w] 
step2  =  D[stepl,  w] 

step3  =  step2/ . kr [w] ->kr/ . (kr" \ [Prime] ) [w] ->dkdw 
ddkdndw  =  Simplify [ step3 ]  (*  d~2  kr  /  (  dn  dw)  *) 

beta  =  — ( 1 / w)  dkdn/ddkdndw  /. (kr/w)  ->Sp 
Simplify [%] 

Equation  (4.26)  is  valid  when  the  bottom  is  f  uid  halfspace,  as  is  the  case  in  all  of 
the  examples  in  this  chapter.  One  could  extend  the  present  analysis  to  include  the 
effects  of  an  arbitrarily  layered  (but  lossless)  bottom  by  letting  the  bottom  ref  ection 
coeffcient  be  a  function  of  both  kr  and  uj  independently:  ol/Ay,  cr).  This  can  be 
done  by  changing  the  f  rst  line  of  the  above  Mathematica  code  to: 

ff  =  n  -  (w  phil [kr/w]  +  phidn[kr,w]  +  1) 
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The  resulting  expression  is  extremely  complicated,  but  analytic  none-the-less. 
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Chapter  5 

Array  processing  with  the 
waveguide  invariant 


Chapter  3  used  striations  in  I  (r.  uj  )  measured  from  a  single  hydrophone  to  estimate 
the  range  to  an  acoustic  source,  by  assuming  that  (5  =  1.  Chapter  4  showed  that 
surface-ref  ected  bottom-ref  ected  (SRBR  —  see  Sec.  2.2.2)  modes  usually  had  (3mi 
values  close  to  one,  whereas  non-SRBR  modes  had  (3mn  values  that  varied  widely 
and  depended  strongly  on  the  details  of  the  SSR  In  this  chapter,  we  develop  a 
technique  for  observing  waveguide  invariant  striations  with  a  horizontal  line  array 
(HLA),  which  allows  for  rejecting  noise  from  sources  other  than  the  source  of  in¬ 
terest.  We  also  show  that  an  HLA  with  suff  cient  aperture  can  f  Iter  out  non-SRBR 
modes,  resulting  in  a  striation  pattern  consisting  only  of  SRBR  modes  which  are 
likely  to  have  /3mi  «  1. 

As  with  the  rest  of  this  thesis,  mode  numbers  are  denoted  by  m  and  l. 


5.1  Introduction 

Passive  source  localization  is  one  of  the  most  straightforward  applications  of  the 
waveguide  invariant.  If  one  uses  a  single  hydrophone  to  measure  the  acoustic  inten¬ 
sity  in  a  shallow-water  waveguide  while  a  broadband  source  moves  away  from  the 
hydrophone,  the  resulting  plot  of  acoustic  intensity  will  contain  striations  described 
by  the  waveguide  invariant.  The  slopes  of  the  striations  can  be  used  to  estimate  the 
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range  to  the  source,  as  was  done  in  Ch.  3  and  in  [56],  This  type  of  processing  is  re¬ 
ferred  to  as  incoherent  processing  because  it  depends  only  on  the  acoustic  intensity 
(as  opposed  to  the  complex  pressure  level)  and  is  discussed  further  in  Sec.  5.2. 

Although  observing  intensity  striations  using  a  single  hydrophone  is  appealing 
because  of  its  simplicity,  the  disadvantage  is  that  if  there  are  multiple  sources  in  the 
water,  all  of  those  sources  contribute  to  the  striation  pattern  in  a  non-linear  manner 
and  the  striations  from  the  source  of  interest  can  become  obscured.  To  combat 
this  problem,  one  can  use  an  array  and  beamform  towards  the  source  of  interest. 
The  use  of  an  array  to  observe  the  striations  necessitates  coherent  processing  of  the 
measured  acoustic  f  eld. 

In  [66],  Yang  showed  that  the  output  of  a  horizontal  line  array  (HLA)  with  a 
uniform  array  weights  can  be  used  to  observe  waveguide  invariant  striations.  The 
output  of  the  beamformer  will  contain  striations  that  are  almost  the  same  as  the 
striations  present  in  the  single-hydrophone  acoustic  intensity,  while  rejecting  noise 
coming  from  directions  other  than  that  of  the  source  of  interest.  In  [62],  Turgut  et 
al.  used  the  methods  described  in  [66]  to  successfully  localize  an  acoustic  source 
with  experimental  data. 

In  Sec.  5.3,  it  is  shown  that  when  the  HLA  is  short  and  the  temporal  frequency 
is  low  (small  aperture  relative  to  wavelength),  processing  the  array  data  with  a 
uniformly  weighted  beamformer  will  yield  the  expected  results:  the  array  processor 
output  will  contain  desired  striations  from  the  source  of  interest  while  rejecting 
noise  from  other  directions.  But  when  the  array  is  long  or  the  temporal  frequency 
is  high  (large  aperture  relative  to  wavelength),  uniformly  weighted  beamforming 
can  have  an  undesired  effect:  the  main  lobe  of  a  uniformly  weighted  beamformer 
can  get  so  narrow  that  it  eliminates  some  of  the  modes  necessary  to  create  the  modal 
interference  pattern  (striations)  that  is  one  is  trying  to  observe. 

This  chapter  shows  how  to  design  HLA  array  weights  that  are  optimal  for  ob¬ 
serving  waveguide  invariant  striations  from  a  single  source  while  rejecting  noise 
from  other  directions.  This  chapter  also  analyzes  when  it  is  benef  cial  to  design 
such  array  weights  as  opposed  to  simply  using  uniform  ones. 

Section  5.4  uses  simulated  data  to  illustrate  how  an  HLA  can  f  Iter  out  non- 
SRBR  modes  which  typically  don’t  have  (3mi  «  1.  Sec.  5.5  uses  experimental 
data  (from  the  same  experiment  that  was  used  in  Ch.  3)  to  show  that  an  HLA  can 
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reject  noise  but  preserve  the  waveguide  invariant  striations.  Sec.  5.6  summarizes 
and  concludes  the  chapter. 

Before  moving  forward  with  the  present  analysis,  it  is  worth  noting  that  other 
array  processing  techniques  have  been  developed  to  exploit  the  waveguide  invari¬ 
ant,  including  source  localization  using  broadband  matched  f  eld  processing  [59], 
and  covariance  matrix  estimation  for  adaptive  beamforming  [57,  37]  (See  Appendix 
A).  In  the  present  analysis,  the  focus  is  on  using  an  HLA  to  f  Iter  out  non-SRBR 
modes,  and/or  observe  the  waveguide  invariant  striations  from  one  source  while 
eliminating  contributions  to  the  acoustic  intensity  from  other  sources. 

5.2  Review  of  incoherent  processing 

As  discussed  in  Sec.  2.3.1,  the  acoustic  intensity  in  a  waveguide  resulting  from 
a  single,  compact  acoustic  source  can  be  written  as  summation  of  cosines  [see 
Eq.  (2.10)].  Many  processing  schemes  that  utilize  the  waveguide  invariant,  such 
as  the  one  presented  in  Ch.  3,  rely  on  accurately  observing  the  striations  resulting 
from  the  cosine  terms  in  Eq.  (2.10).  But  when  other  sound  sources  are  present 
in  the  ocean,  the  striations  from  the  source  of  interest  can  be  become  diff  cult  to 
observe  amongst  all  the  noise.  An  example  of  this  can  be  seen  in  Fig.  5.4(a). 

If  one  has  an  array  of  hydrophones  rather  than  a  single  hydrophone,  the  stri¬ 
ations  from  the  source  of  interest  can  be  enhanced  by  using  coherent  processing. 
But  doing  so  requires  a  judicious  choice  of  array  weights,  as  will  be  discussed  in 
the  next  section. 


5.3  Coherent  processing 

The  concept  of  observing  waveguide  invariant  striations  using  an  HLA  was  origi¬ 
nally  introduced  by  Yang  in  [66],  which  analyzed  the  special  case  of  an  HLA  with 
uniform  array  weights.  In  this  section  we  generalize  Yang’s  results  to  array  pro¬ 
cessors  (beamformers)  with  non-uniform  weights,  such  as  a  Hamming  Window  or 
weights  designed  by  more  advanced  techniques.  We  restrict  our  attention  to  uni¬ 
formly  spaced  HLAs.  It  will  be  shown  that  the  trade-off  to  be  considered  when 
choosing  the  weights  is  not  the  usual  “beam  width  versus  side  lobe  level”  that  one 
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encounters  with  planewave  array  processing,  because  the  desired  signal  is  not  a 
planewave  but  rather  the  striations  resulting  from  the  interference  between  a  set  of 
modes. 

We  start  by  def  ning  some  notation  (bold  faced  lower-case  letters  denote  vector 
quantities): 

•  co  is  the  speed  of  sound  at  the  HLA  depth. 

•  ko  —  ~  is  the  wavenumber  at  the  HLA  depth. 

•  J  is  the  number  of  elements  in  the  HLA,  with  each  element  numbered  as 
3  =  1 ,  3  =  2,  j  =  J. 

•  is  the  look  direction  as  measured  from  broadside  of  the  HLA. 

•  kL  =  A'o  sin  tfr  is  the  horizontal  wavenumber  corresponding  to  the  look  di¬ 
rection. 

•  v(A;l)  =  [ui  (&£,),  U2  (&£,),  •••  ,  vj(kL)]T  is  the  steering  vector  correspond¬ 
ing  to  the  look  direction,  v  is  a  column  vector. 

•  w  =  [w\,  W2,  ■  ■  ■  ,  wj]T  the  array  weights  column  vector  (e.g.,  Hamming, 
Hann,  uniform),  w  is  a  real-valued  vector  and  does  not  include  the  phase 
shift  of  the  steering  vector. 

•  9s  is  the  bearing  of  the  acoustic  source  as  measured  from  broadside  of  the 
HLA,  in  the  plane  parallel  to  the  ocean  surface. 

•  T  is  the  wavenumber  response  function  of  the  array  processor,  which  de¬ 
pends  on  the  array  geometry  and  the  array  weights  w,  and  is  described  in 
more  detail  later. 

5.3.1  HLA  response  to  a  planewave 

An  array  and  its  weights  can  be  characterized  by  their  wavenumber  response  func¬ 
tion,  denoted  here  by  T  [61,  Chs.  2  and  3],  T  describes  the  gain  of  an  array  proces¬ 
sor  for  a  planewave  with  a  given  wavenumber  along  the  axis  of  the  array.  As  shown 
in  [61,  Eq.  2.122],  T  for  a  uniformly  spaced  HLA  depends  only  on  the  difference 
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between  wavenumber  corresponding  to  look  direction  (, )  and  the  wavenumber 
corresponding  to  the  source  direction  (ks),  as  opposed  to  depending  on  hr,  and  ks 
independently.  T  also  depends  on  (is  parameterized  by)  the  array  weights.  Thus 
we  write  the  wavenumber  response  function  as  T(kL  —  ks',  w).  It  is  well  known 
that  when  w  =  1,  T {Ul  —  ks ;  1)  is  a  periodic  sine  function  of  {hi,  —  ks)  [61, 
Eq.  2.94]. 

Let  a  superscript  H  represent  a  complex  conjugate  transpose,  and  an  asterisk 
represent  the  complex  conjugate  of  a  scalar.  The  wavenumber  response  function  is 
then  def  ned  as  [61,  Eq.  2.58] 

J 

T (kL  -  ks;  w)  =  v(ks)H(v(kL) w)  =  'Y^vj(ks)*(vj(kL)wj)  (5.1) 

3= 1 

where  (v(fc^)w)  is  the  element-by-element  multiplication  of  v(/,yj  and  w. 

For  an  HLA  with  J  elements  numbered  1  to  J  and  an  inter-element  spacing  of 
d,  the  j-th  element  is  located  at  a  distance  of  d(j  —  1)  from  the  1st  element,  and 
the  center  of  the  array  is  at  a  distance  of  d{J  —  l)/2  from  the  1st  element.  Let 
dj  represent  the  displacement  of  the  j-th  element  of  the  array  to  the  center  of  the 
array: 

dj=d(j  -l)-d(J-l)/2  (5.2) 

The  j-th  element  of  the  steering  vector  v  is  then  (Eq.  2.78  of  Van  Trees) 

Vj(kL)  =  exp  [+ ikL,dj\  (5.3) 


In  this  chapter  we  only  consider  w  that  are  real  and  symmetric  with  respect  to 
the  center  of  the  array.  Because  of  that,  and  because  dj  was  def  ned  with  respect  to 
the  center  of  the  array,  T (kL  —  ks;  w)  will  be  real.  This  is  a  result  of  the  Discrete 
Fourier  Transform  relationship  between  T  and  w  [61,  Sec.  2.3].  Knowing  that 
T  ( k  fJ  —  ks;  w)  is  real  will  help  simplify  the  algebra  later  in  this  chapter,  but  is  of 
little  conceptual  importance. 
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5.3.2  HLA  response  to  source  in  waveguide 

T(kL— ks;  w)  represents  the  response  of  an  array  to  the  pressure  f  eld  of  a  planewave. 
In  this  subsection  we  derive  the  response  of  an  HLA  array  with  weights  w  to  the 
pressure  f  eld  from  a  discrete  source  in  a  waveguide.  The  intensity  of  that  response, 
when  plotted  as  a  function  of  range  and  frequency,  will  contain  waveguide  invariant 
striations.  [66]  contains  a  similar  derivation  for  the  specif  c  case  of  uniform  array 
weights  (w  =  1).  In  the  present  work,  we  generalize  Yang’s  result  and  explain 
some  advantages  of  using  non-uniform  array  weights. 

Let  T total (w j  0l,0s,  rc:  w)  represent  the  response  function  of  an  HLA  due  to  a 
source  at  a  bearing  of  Os  and  a  range  of  rc  in  a  range  independent  waveguide,  while 
looking  in  a  direction  Ol-  In  contrast  to  T  which  depends  only  on  ks  —  ks,  Ttotai 
depends  on  uj,  Ol  and  Os  for  reasons  that  will  become  clear  later;  it  also  depends 
on  the  array  weights  w.  It  will  be  shown  that  Ttotai  can  be  written  as  a  sum  of  Ts, 
one  for  each  mode. 


Ttotaifw,  0L ,  0S,  rc:  w)  =  pH(v(kL)w)  (5.4) 

where  p  is  the  pressure  at  the  array  elements  due  to  the  source  in  the  waveguide. 
Note  that  ks  is  a  function  of  Ol,  as  def  ned  at  the  beginning  of  this  section. 

To  calculate  p,  let  rc  be  the  range  from  the  source  to  the  center  of  the  HLA.  If 
the  distance  between  the  array  and  the  source  is  large  compare  to  the  array’s  length, 
then  the  range  from  the  source  to  the  j- th  element  of  the  array  can  be  approximated 
as 


rc  +  sm0sdj  (5.5) 

The  pressure  at  the  j-th  element  of  the  array  can  written  by  replacing  the  r  in 
Eq.  2.7  with  Eq.  5.5,  resulting  in: 

Pi  =  E  Am  exp  [+ik  rmix c  +  sin  Osdj)\.  (5.6) 

m 

Eq.  (5.6)  can  be  interpreted  as  a  sum  of  planewaves  along  the  array,  each  with 


114 


amplitude  an  Am,  a  wavenumber  along  the  array  of  krm  sin  Os,  and  a  phase  of 

Re¬ 


inserting  Eqs.  (5.6)  and  (5.3)  into  Eq.  (5.4),  and  rearranging  the  summations 
yields: 


Ttotai(w;0L,0s,Rc;  w)  =  p'Hv(fcL)w 
J 


(5.7) 


^2PjVj(kL)wj 


*j\  wj 


3= 1 

^  I  ^  Am  exp  [+ikrm(rc  +  sin  9sdj)}  j  exp  [+ikLdj 
j= l  \  m  J 

J 

exp  [— i&rmrc]  ^  exp  [— ikrm  sin  9sdj]  exp  [+ikidj]wj  (5.8) 

m  j= i 


The  summation  over  j  in  Eq.  (5.8)  is  the  same  as  Eq.  (5.1),  but  with  the  ks  in 
Eq.  (5.1)  replaced  by  krm  sin  9s-  Eq.  (5.1)  is  the  defnition  of  T  for  a  planewave. 
We  can  therefore  write  Ttotai  as  a  sum  of  Ts  -  one  for  each  mode  m. 


Def  ne  kms  =  krm  sin  6s,  which  is  the  wavenumber  of  mode  m  along  the  array 
due  to  a  source  at  a  bearing  9s-  Ttotai  can  then  be  written  as 

Ttotai (w,  Ol,  Os',  w)  =  Am  exp  [-ikrmrc]T(uj,  kL  -  kmS)  (5.9) 

m 

bearing  in  mind  that  kL  is  a  function  of  Ol,  and  kms  is  a  function  of  9s-  To  the 
HLA,  each  mode  appears  to  be  a  planewave  propagating  along  the  array  with  a 
wavenumber  of  kms,  with  an  amplitude  of  Am,  and  a  phase  shift  of  krmrc.  So 
Ttotai  is  a  sum  of  Ts  corresponding  to  each  mode. 


The  magnitude  of  Ttotai(w,  Ol,  Os',  w)  contains  the  waveguide  invariant  stria- 
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tions  analogous  to  those  in  Eq.  (2.10): 


|Ttotai(w,0L,0s;  w)|2  =  Am  exp  [-ikrmrc]T (u ,  kL  -  kmS; 

m 

=  J2A2q\r(kL-kqS-,™)  |2 


W, 


Q 

+  ^  ArnAicos(Akmirc)(T(kL  -  kmS]  w))(Y(kL  -  kiS;  w)) 

ra,Z;  m^l 


(5.10) 


where  A kmi  =  krrn  —  kri  is  the  difference  in  the  horizontal  wavenumbers  between 
mode  m  and  mode  l,  and  we  have  used  the  fact  that  T (kL  ~  kms;  w)  is  real. 

A  careful  term-by-term  comparison  of  Eq.  (5.10)  to  Eq.  (2.10)  is  essential  to 
understanding  the  present  analysis.  The  waveguide  invariant  striations  are  a  re¬ 
sult  of  the  cos(A kmirc)  terms,  which  represent  the  interference  between  modes. 
The  cosine  terms  in  Eq.  (5.10)  are  the  same  as  in  Eq.  (2.10)  but  are  multiplied  by 
(T (kL  -  kmS-,w))(Y(kL  -  kiS;  w)). 

In  order  to  “extract”  the  waveguide  invariant  striations  from  HLA  measure¬ 
ments  of  the  acoustic  f  eld,  one  would  f  rst  compute  the  complex  pressure  level  at 
each  hydrophone  of  the  array  at  each  snapshot.  The  resulting  complex  pressure 
levels  would  then  be  fed  into  an  array  processor,  and  the  magnitude  of  the  output 
of  the  array  processor  would  be  plotted  versus  frequency  and  snapshot  number. 

For  the  purpose  of  understanding  the  present  work,  it  is  easiest  to  consider 
a  scenario  where  the  source  changes  range  but  not  bearing  between  snapshots  so 
that  one  can  plot  |Ytotai(u;,  0l,  Os',  w)|2  versus  rc  and  u  while  keeping  0l  f  xed  at 
0s-  (In  a  practical  situation,  one  would  likely  plot  |Ttotai(w,  0l,  0s]  w)|2  versus 
snapshot  number  and  u  while  changing  0  jJ  with  the  snapshot  number  to  follow  the 
source,  and  then  estimate  the  range  and  speed  of  the  source.  See  [62].) 


5.3.3  Array  processing  for  the  waveguide  invariant 

To  simplify  notation,  let  A k  equal  the  difference  between  a  wavenumber  received 
by  an  array,  and  the  wavenumber  corresponding  to  the  look  direction  of  array  pro¬ 
cessor.  That  is,  A k  =  kL  —  ks  in  free  space  or  A k  =  kL  —  km.s  in  a  waveguide. 
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The  purpose  of  array  processing  is  to  exclude  contributions  to  the  output  of  the 
array  processor  from  sources  located  at  a  bearings  different  than  that  of  the  source 
of  interest.  Traditionally,  array  processing  (beamforming)  was  developed  to  distin¬ 
guish  planewaves  from  one  another,  and  thus  the  goal  was  the  design  array  weights 
that  resulted  in  as  narrow  of  a  main  lobe  in  T  as  is  possible  while  controlling  the 
level  of  the  side  lobes  in  T  [61,  Ch.  3].  The  array  weights  were  usually  chosen  so 
that  Y(A k;  w)  was  largest  at  Ak  =  0  and  T(A k;  w)  generally  decreased  in  value 
as  A k  increased. 

In  contrast  to  the  planewave  model,  the  waveguide  invariant  striations  are  the 
result  of  modes  interfering  with  one  another.  In  the  previous  subsection  it  was 
shown  that  each  mode  appears  to  the  HLA  as  a  planewave  with  wavenumber  kms ■ 
If  the  array  processor  weights  have  a  high  enough  wavenumber  resolution  (i.e., 
main  lobe  width)  to  distinguish  a  low  order  mode  (e.g.,  kis)  from  a  high  or¬ 
der  model  (e.g.,  kgs)  coming  from  the  same  source,  then  the  striations  result¬ 
ing  from  the  interference  term  between  the  two  modes,  cos  ( A kmi  rc)(T(k l  — 
kmS'i  w))(T(fc^/  —  kis;  w))  from  Eq.  (5.10),  will  have  a  very  small  magnitude  be¬ 
cause  either  (T(/j£,  —  kms;  w))  or  (T (&l  —  kis', w))  will  be  small  no  matter  what 
the  value  of  ks  is.  That  is,  no  matter  where  the  array  is  steered,  the  desired  stria¬ 
tions  term  may  be  suppressed  due  to  the  high  resolution  of  the  array  processor.  We 
investigate  this  further  in  the  next  subsection,  and  then  suggest  how  to  design  better 
array  weights  later  in  this  chapter. 


5.3.4  Uniform  weighs 


This  subsection  analyzes  the  circumstances  which  result  in  suppressed  striations  for 
the  case  of  an  array  processor  with  uniform  weights.  [66]  derived  T total  (bb  0l,9s',  w) 
for  the  specif  c  case  of  w  =  1.  In  that  case,  T  is  a  periodic  sine  function  [61, 
Eq.  2.94]: 


T(A k]  w) 


sin  (4pA k) 
sin  (| Ak) 


(5.11) 
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Inserting  this  into  Eq.(5.10)  yields 


<? 


sin  (4r(fcL  -  kqs)) 
sin  -  A:gs)) 


+  ^2  AmAi  cos(A fcmZrc) 

ra,Z;  m^l 


(  sin  (^(^l  ~  fcmg))\ 
\  sin  -  kmS))  ) 


(5.12) 

sin  (^(^l  ~  hs))\ 
sin  {^(kL  —  kis))  ) 


Eq.  (5.12)  is  equivalent  to  [66,  Eqs.  14  and  15].  The  second  term  of  Eq.  (5.12) 
contains  the  same  cosine  terms  that  are  in  Eq.  (2.10),  but  they  are  multiplied  by 
two  periodic  sine  functions. 

When  (\Jd/2Ak\  <C  7 r),  the  periodic  sine  function  is  near  its  maximum  value. 
As  |  Jd/2Ak\  increases  in  value,  the  periodic  sine  function  decreases  in  value  until 
|  Jd/2Ak\  >  7 r,  at  which  point  the  periodic  sine  function  will  have  a  very  small 
value  (i.e.,  be  in  the  side  lobe  of  the  periodic  sine). 

(Jd/2)  is  approximately  half  of  the  array  length.  So  if  the  array  length  times 
( kis—kms )  differ  by  an  amount  larger  than  twice  the  main  lobe  width  of  T(Afc;  w), 
then  the  cosine  term  in  Eq.  (5.12)  is  going  to  be  very  small  no  matter  what  the  value 
of  kL  is. 

This  can  occur  in  practical  situations.  Consider  a  waveguide  where  the  water 
column  has  a  sound  speed  of  1500  m/s  and  the  seabed  has  a  sound  speed  of  1600 
m/s.  If  a  source  with  a  center  frequency  of  500  Hz  is  at  end  f  re,  then:  kms  =  krm, 
the  frst  order  mode  has  kr\  «  27r500/1500,  and  highest  order  mode  has  krM  ~ 
27r500/1600.  An  array  that  is  longer  than  about  200  meters  (with  w  =  1  )  will 
have  such  a  narrow  main  lobe  that  there  is  no  choice  of  k  jJ  for  which  T  for  mode  1 


sin  (4p  (kL  ~  kr i)) 
sin  (f  (kL  -  kr i)) 


and  T  for  mode  M 

sin  (f  (kL  -  krM))  (5  14) 

sin  (|  (kL  ~  krM )) 
are  in  their  respective  main  lobes. 

Although  one  may  not  always  want  the  term  corresponding  to  Ak\M  to  be 
in  the  array  output,  an  analogous  argument  can  be  made  for  any  two  modes  with 
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differing  horizontal  wavenumbers.  Even  if  the  two  modes  are  both  in  the  main 
lobe,  their  cosine  term  will  be  suppressed  to  some  extent  because  the  periodic  sine 
function  starts  to  decrease  in  value  as  soon  as  the  wavenumber  difference  in  the 
argument  isn’t  zero. 

In  this  subsection,  we  showed  why  uniform  weights  (w  =  1  )  can  have  unde¬ 
sired  results.  In  the  next  subsection,  we  show  how  to  choose  the  weights  to  prevent 
the  suppression  of  the  desired  cosine  terms. 

In  addition  to  the  issue  analyzed  in  this  section  regarding  the  suppression  of 
cosine  terms  in  the  summation,  there  is  one  other  potential  problem  with  the  array 
processor  output.  As  noted  in  [66],  ( k,L  —  kms )  and  ( k,L  —  hs )  are  functions  of  c o. 
And  so  multiplying  the  cosine  term  with 

sin  (4t(A;l  -  kmS))  sin  (ff  (fcL  -  kiS)) 
sin  (| (kL  -  kmS))  sin  (f  (kL  -  kiS )) 

may  have  an  adverse  affect  on  the  striation  pattern,  which  lies  in  the  (rc,  u)  plane. 
It  will  become  clear  in  the  next  subsection  that  designing  a  w  that  prevents  the 
suppression  of  the  cosine  terms  will  also  prevent  T{kjJ  —  kms-  w)  and  T(k^  — 
kis',  w)  from  being  a  strong  function  of  frequency  and  thus  preserve  the  desired 
striations  resulting  from  the  cosine  terms. 

5.3.5  Designing  array  weights  for  waveguide  invariant  striations 

The  previous  subsection  discussed  why  uniform  weights  can  lead  to  a  suppression 
of  some  striation  terms  because  uniform  weights  correspond  to  a  narrow  main  lobe 
in  T.  In  this  subsection  we  discuss  how  to  choose  non-uniform  weights  for  the 
purpose  of  observing  the  waveguide  invariant  striations.  To  understand  this  subsec¬ 
tion,  it  is  useful  to  think  of  the  array  processor  as  a  spatial  f  Iter  whose  wavenumber 
response  is  H(k )  where  k  is  the  wavenumber  along  the  array. 

In  the  case  of  planewave  beamforming  —  where  the  goal  is  to  estimate  the 
signal  coming  from  a  one  direction  while  rejecting  signal  coming  from  other  direc¬ 
tions  —  the  ideal  H(k)  is  an  impulse  in  k  space.  Because  that  cannot  be  achieved 
with  a  f  nite  aperture  array,  one  chooses  a  trade-off  main  lobe  width  and  side  lobe 
levels. 
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But  in  the  case  of  using  an  array  to  observe  the  waveguide  invariant  striations, 
the  ideal  H  (. k )  is  not  an  impulse  in  k  space;  the  interference  pattern  is  the  result 
of  modes  with  different  k  values  interfering  with  each  other,  and  so  a  range  of  k 
values  must  pass  through  the  array  processor. 

An  ideal  H (k)  for  observing  the  waveguide  invariant  striations  is  an  H(k ) 
that  passes  through  all  wavenumbers  associated  with  the  cosine  terms  in  Eq.(5.10) 
while  blocking  out  all  other  wavenumbers,  so  that  none  of  the  striation  terms  are 
suppressed  but  noise  coming  from  directions  other  than  the  source  direction  get 
blocked  out  (other  ideal  H (k) s  are  discussed  later  in  this  subsection).  For  example, 
in  a  Pekeris  waveguide  with  cwater  =  1500  m/s  and  Cbottom  =  1600  m/s  H(k),  all 
modal  horizontal  wavenumbers  will  be  between  y^  and  y^ .  So  the  array  proces¬ 
sor  should  pass  all  wavenumbers  between  k  =  sin  Os  ( y^j )  and  k  =  sin  Qg  (y^ ) 
with  equal  amplitudes  while  blocking  signals  with  all  other  wavenumbers: 


H ideal  {k) 


0  k  <C  sin  Os  (  1(500  )  =  ^lower 

<  1  sin °S  (lio)  <  k  < sin °s  (tsso) 
0  k  >  sin  0s  (  1500  )  =  Supper 


(5.16) 


//ideal  (k)  in  Eq.  (5.16)  is  a  complex  band-pass  spatial  flter  whose  spatial  cutoff 
wavenumbers  depend  on  the  temporal  frequency  oj  and  the  source  direction  9$.  It 
is  complex  in  the  sense  that  the  input  signal  is  complex  (the  complex  pressure  level 
along  the  array  at  some  value  of  u).  So  the  flter  must  treat  positive  and  nega¬ 
tive  wavenumbers  differently,  in  contrast  to  typical  time  domain  band-pass  f  Iters 
whose  impulse  response  is  generally  conjugate  symmetric  (which  treat  positive  and 
negative  frequencies  the  same  when  the  input  signal  is  real). 

H ideal  (k)  can  only  be  realized  with  a  inf  nitely  long  array,  so  it  is  necessary  to 
approximate  H^eai(k).  This  can  be  done  with  a  f  nite  impulse  response  (FIR)  flter 
of  length  J  because  there  are  J  elements  in  the  array.  A  straightforward  approach  to 
designing  the  FIR  flter  is  to  f  rst  design  a  linear-phase,  low-pass  FIR  flter  centered 
at  k  =  0  using  standard  digital  flter  design  techniques  [46,  Ch.  7].  The  resulting 
FIR  flter  will  be  the  array  weights  w.  (Note  that  a  linear  phase  flter  is  required 
because  it  yields  a  symmetric  w).  Then  multiply  w  by  the  steering  vector  v  to  shift 
the  pass-band  from  being  centered  at  k  =  0  to  being  centered  at  some  k  value  that 
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is  halfway  between  the  fciower  and  /cupper  in  Eq.  (5.16). 

For  example,  in  the  above  case,  one  would  f  rst  design  a  low-pass  f  Iter  with  a 
cutoff  frequency  of  fccutoff  =  (Supper  —  Mower )/2.  One  would  then  multiply  that  f  Iter 
by  v(fccenter)  to  shift  the  center  frequency  from  up  to  kcen ter  =  (kupper  +  howeT)/2. 
The  resulting  output  can  then  be  described  by  Eq.  (5.10)  where  the  &L  =  Center- 
Designing  w  involves  choosing  between  many  trade-offs,  which  are  discussed  in 
digital  f  Iter  design  text  books  such  as  [46],  The  next  section  discusses  one  possible 
choice  of  w  in  depth. 

iTideai(fc)  depends  on  the  Os  and  the  sound  speeds  in  the  environment,  which 
may  not  be  known  accurately.  In  that  case,  the  pass-band  of  the  flter  can  be 
widened  at  the  expense  of  letting  in  more  noise.  Or  one  could  narrow  the  pass- 
band  in  order  to  reject  more  noise,  at  the  expense  of  suppressing  some  of  the  cosine 
terms. 

In  some  cases,  one  may  purposely  choose  a  flter  that  is  narrower  than  the 
-^ideai(^)  in  Eq.  (5.16).  The  waveguide  invariant  striations  (cosine  terms)  are  most 
useful  when  the  value  of  8mi  is  predictable  (i.e.,  not  sensitive  to  the  details  of  the 
environment).  But  8ml  is  only  predictable  under  certain  circumstances:  Sec.  2.3 
and  Ch.  4  qualitatively  showed  that  when  both  modes  are  SRBR  modes  that  are 
not  near  cutoff,  8mi  «  1.  For  other  modes,  the  value  of  8  is  often  unpredictable. 
Therefore,  it  may  be  useful  to  flter  out  wavenumbers  that  correspond  to  modes 
close  to  cutoff,  and  to  flter  out  wavenumbers  that  correspond  to  non-SRBR  modes. 
Then  the  cos  (A kmir)  terms  with  unpredictable  8ml  values  will  be  suppressed  and 
the  array  processor  output  may  be  more  useful. 

It  is  worth  noting  that  f  ltering  out  wavenumbers  corresponding  to  non-SRBR 
modes  while  keeping  those  corresponding  to  SRBR  modes  requires  fairly  high 
wavenumber  resolution,  and  thus  requires  a  long  array  relative  to  the  wavenum¬ 
ber  differences. 


5.4  A  simulated  example 

To  validate  the  theory  presented  in  the  previous  section  and  explain  some  of  the 
details  required  for  implementation,  a  simulated  example  is  provided  here.  The 
environment  and  array  parameters  are  chosen  specif  cally  to  illustrate  the  concepts 
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from  the  previous  section,  but  are  realistic  none-the-less. 

The  environment  is  a  range-independent  shallow-water  waveguide  with  a  typi¬ 
cal  summer  sound  speed  prof  le  (SSP),  shown  in  Fig.  5. 1.  The  bottom  is  a  halfspace 
with  sound  speed  of  1600  m/s,  a  density  of  1 .5  g/cm3,  and  a  attenuation  of  1/2  dB/A. 
The  array  was  a  256  element  F1LA  with  a  spacing  corresponding  to  A/2  at  700  Hz; 
the  total  array  length  was  273  meters.  The  source  was  at  endf  re,  which  maximizes 
the  effective  aperture  of  the  array. 

The  acoustic  f  eld  was  calculated  using  Kraken  [48]  from  a  range  of  3000  m  to 
5000  m,  and  a  temporal  frequency  of  300  Hz  to  700  Hz.  The  source  and  receiver 
depth  (52  and  84  meters,  respectively)  were  chosen  to  be  below  the  thermocline  so 
that  the  bottom-trapped  modes  would  contribute  to  the  total  pressure. 

5.4.1  Filter  design 

As  discussed  in  Sec.  5.3.5,  it  can  be  useful  to  f  Iter  out  modes  that  do  not  interact 
with  the  top  and  bottom  of  the  waveguide,  because  those  modes  are  unlikely  not 
lead  to  a  [3mi  value  that  is  close  to  one.  Modes  that  are  close  to  cutoff  will  also 
not  have  fimi  ~  1,  but  modes  that  are  close  to  cutoff  have  very  small  amplitudes 
because  they  are  highly  attenuated,  so  it  is  usually  not  necessary  to  f  Iter  them  out. 
Therefore  the  ideal  array  processor  for  this  environment  would  f  Iter  all  wavenum¬ 
bers  except  those  with  phase  speeds  between  1535  m/s  and  1600  m/s,  because 
modes  with  those  phase  speeds  are  top  and  bottom  interacting  modes. 

To  create  an  linear-phase  FIR  f  Iter  that  approximates  the  ideal  f  Iter,  a  tech¬ 
nique  that  minimizes  the  integrated  squared  error  between  the  ideal  f  Iter  and  the 
magnitude  response  of  the  flter  was  used  (specif  cally,  the  firls  tool  in  MAT- 
LAB).  Because  the  wavenumber  associated  with  a  specif  c  phase  speed  is  a  func¬ 
tion  of  frequency  (e.g.,  k  =  j/j/y),  a  different  FIR  flter  must  be  designed  for  each 
temporal  frequency.  The  width  of  the  pass-band  is  A k  =  u( so  at 
higher  frequencies  the  FIR  flter  must  have  a  wider  pass-band. 

5.4.2  Simulated  results 

Figure  5.2  shows  the  response  of  the  FIR  f  Iters  (array  processor  weights)  designed 
for  (a)  /  =  300  Hz,  (b)  /  =  500  Hz,  and  (c)  and  /  =  700  Hz.  Also  shown  is  the 
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response  of  uniform  array  weights,  which  is  the  same  for  all  temporal  frequencies. 
In  subf  gure  (a)  (/  =  300  Hz),  the  FIR  f  Iter  response  is  similar  to  the  uniform 
weights,  but  slightly  wider.  In  subf  gure  (c)  (/  =  700  Hz),  the  difference  between 
the  FIR  f  Iter  and  the  uniform  weights  is  very  apparent;  the  FIR  f  Iter  has  a  wide, 
f  at,  pass-band  that  allows  all  the  desired  modes  the  pass  through. 

As  discussed  in  Sec.  2.3,  waveguide  invariant  striations  are  of  the  form 


SJ_ 

5r 


(5.17) 


where  ^  is  the  slope  of  the  striation  and  ft  is  the  waveguide  invariant.  So  striations 
with  /3  =  1  correspond  to  lines  of  the  form  /  =  Cft  r  where  (ft  is  an  arbitrary 
constant.  Figure  5.3  shows  the  intensity  from  the  (a)  unf  ltered  pressure  and  (b) 
array  processed  (f  ltered)  pressure.  The  black  lines  correspond  to  (3  =  1.  Notice 
that  the  striations  in  subf  gure  (b)  correspond  much  more  closely  io  ft  =  1  than 
do  the  striations  in  subf  gure  (a),  because  the  pressure  from  which  the  intensity  in 
subf  gure  (b)  was  calculated  was  f  ltered  (array  processed)  to  remove  modes  that 
were  not  top  and  bottom  interacting. 


SSP 


Figure  5.1:  The  sound  speed  prof  le  for  the  simulated  environment.  It  is  a  simplif  ed 
model  of  a  shallow-water  waveguide  in  the  summer.  The  bottom  halfspace  starts  at 
100  meters. 
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(a)  (f  =  300  Hz) 


(b)  (f  =  500  Hz) 


(c)  (f  =  700  Hz) 


Figure  5.2:  Plots  of  the  wavenumber  response  function  of  three  different  array 
weights,  each  designed  for  a  different  temporal  frequency  /.  The  response  due 
to  uniform  array  weights  is  plotted  in  each  subf  gure  as  well,  for  comparison,  (a) 
/  =  300  Flz.  (b)  /  =  500  Flz.  (c)  /  =  700  Flz.  As  /  increases,  the  array  weights 
change  so  that  the  main  lobe  width  increases  to  ensure  that  all  the  desired  modes 
pass  through  the  array  processor  without  being  attenuated.  In  contrast  to  the  array 
weights  designed  specif  cally  for  observing  the  waveguide  invariant,  the  uniform 
array  weights  do  not  change  with  /  and  thus  their  response  function  is  the  same  for 
all  /■ 
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(a)  unfiltered  (b)  filtered 


range  (m)  range  (m) 


Figure  5.3:  Simulated  acoustic  intensity  as  a  function  of  range  and  frequency, 
which  exhibits  striations  described  by  the  waveguide  invariant  /3.  Subfgure  (a) 
is  a  plot  of  the  intensity  from  the  scalar  pressure.  Subf  gure  (b)  is  a  plot  of  the  in¬ 
tensity  from  the  array-processed  pressure  which  had  wavenumbers  corresponding 
to  undesired  modes  f  ltered  out.  The  black  lines  correspond  to  /3  =  1.  The  stria¬ 
tions  subf  gure  (b)  correspond  much  more  closely  to  /3  =  1  than  do  the  striations  in 
subf  gure  (a),  because  the  undesired  modes  have  been  f  ltered  out. 

5.5  Experimental  example 

The  previous  section  showed  an  example  of  how  flter  out  specif  c  wavenumbers 
due  to  a  single  source.  It  required  a  relatively  long  array  in  order  to  have  high 
enough  wavenumber  resolution.  In  this  section,  some  experimental  results  are  pre¬ 
sented  using  a  much  shorter  array.  With  a  short  array,  the  array  weights  designed 
using  the  technique  from  the  previous  section  will  be  very  similar  to  uniform  array 
weights,  because  the  desired  resolution  will  be  close  to  the  highest  resolution  the 
array  is  capable  of,  which  is  achieved  with  uniform  array  weights.  None-the-less, 
the  array  processing  is  very  useful  at  eliminating  sound  coming  from  directions 
other  than  that  of  the  source  of  interest. 

The  data  set  presented  here  is  the  same  experimental  data  that  was  used  for 
Ch.  3.  In  Ch.  3,  only  the  signal  from  a  single  hydrophone  was  processed.  But 
during  the  experiment,  data  was  collected  from  a  horizontal  array  of  hydrophones. 
So  we  now  process  the  same  data  using  the  array  processing  techniques  discussed 
in  this  chapter. 
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An  acoustic  source  was  lowered  to  a  depth  of  40  meters,  and  pseudorandom 
noise  was  projected  from  the  source  over  a  frequency  band  of  350  to  700  Hz.  A  23- 
element,  17.25  m  long  hydrophone  array  was  towed  directly  away  from  the  source 
along  a  line  of  nearly  constant  bathymetry  at  a  speed  of  1.5  m/s.  The  hydrophone 
array  had  an  inter-element  spacing  of  3/4  m,  and  was  towed  at  a  depth  of  50  m.  The 
sea  f  oor  was  about  90  m  deep  with  a  typical  summer  downward  refracting  sound 
speed  prof  le  (SSP).  (See  Sec.  3.4.2  for  more  details). 

Fig.  5.4  (a)  shows  a  plot  of  the  acoustic  intensity  as  measured  by  a  single  hy¬ 
drophone  from  the  array,  and  is  the  exact  same  data  as  that  shown  in  Fig.  3.8. 
Because  13  ~  1  for  the  environment  where  the  data  was  collected  [19],  one  ex¬ 
pects  to  see  waveguide  invariant  striations  from  the  acoustic  source  of  interest  that 
point  toward  the  origin  of  the  coordinate  system  (towards  the  lower  left  of  each 
subfgure).  Although  such  striations  are  clearly  present  at  frequencies  of  350  to 
700  Hz,  there  is  also  noise  present.  One  particularly  loud  noise  source  causes  its 
own  parabolic  shaped  striations  at  all  frequencies.  Figure  5.4  (b)  shows  the  output 
of  the  array  processor,  which  is  using  weights  designed  by  the  method  described  in 
the  previous  section.  Because  the  array  is  short  relative  to  the  desired  wavenumber 
resolution,  the  weights  are  nearly  uniform.  The  striations  from  the  source  of  in¬ 
terest  are  much  more  clearly  visible  in  subf  gure  (b)  than  in  subf  gure  (a),  because 
much  of  the  noise  has  been  f  ltered  out.  Other  array  processing  (not  shown  here) 
conf  rmed  the  presence  of  the  interfering  sound  source  seen  in  subf  gure  (a). 


5.6  Chapter  summary 

Observing  intensity  striations  from  the  waveguide  invariant  is  useful  for  many 
acoustic  sensing  applications,  such  as  the  passive  range  estimation  in  Ch.  3  and 
many  of  the  applications  discussed  in  Appendix  A.  Although  this  can  be  done  with 
a  single  hydrophone,  an  array  can  be  used  to  suppress  noise  coming  from  directions 
other  than  that  of  the  source  of  interest.  In  [66],  Yang  showed  that  the  output  of  an 
array  processor  with  uniform  array  weights  contains  striations  that  are  similar  the 
striations  from  a  single  hydrophone.  The  present  chapter  generalized  Yang’s  result 
for  uniform  weights  to  include  array  processors  with  arbitrary  weights,  and  showed 
why  is  it  benef  cial  to  use  non-uniform  weights. 
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(a)  (b) 
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Figure  5.4:  Acoustic  intensity  from  experimental  data,  plotted  as  a  function  of 
frequency  and  range  from  the  source  of  interest.  If  there  was  no  noise  in  the  data, 
then  all  that  would  be  seen  is  striations  between  350  to  700  Flz  pointing  towards  the 
origin  of  each  subf  gure’s  respective  coordinate  system.  Both  subf  gures  are  from 
the  same  data  set.  Subf  gure  (a)  is  the  acoustic  intensity  as  measured  by  a  single 
hydrophone  from  an  array  towed  away  from  the  source.  The  plot  has  the  expected 
striations  from  the  source  of  interest,  as  well  as  noise  from  other  interfering  sound 
sources  that  were  present  during  the  experiment.  The  parabolic  shaped  striations 
between  1200  and  1800  m  are  from  an  interfering  sound  source.  Subf  gure  (b)  is 
the  same  as  (a),  but  data  from  all  elements  of  the  array  were  used  to  f  Iter  out  noise 
coming  from  directions  other  than  that  of  the  source  of  interest. 
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Because  the  waveguide  invariant  striations  are  the  result  of  modes  interfering 
with  one  another,  designing  array  processing  weights  for  the  purpose  of  observing 
waveguide  invariant  striations  is  different  than  designing  array  processing  weights 
for  planewave  beamforming.  It  was  shown  that  when  an  HLA  is  long  enough 
to  distinguish  the  horizontal  wavenumbers  of  low  order  modes  from  the  horizon¬ 
tal  wavenumbers  high  order  modes,  an  array  processor  with  uniform  weights  can 
have  undesired  effects.  An  example  was  given  demonstrating  how  to  design  array 
weights  that  allow  a  range  of  wavenumbers  to  pass  through  the  array  processor, 
preserving  all  the  desired  striation  terms  while  suppressing  noise. 

It  was  also  shown  that  even  in  the  absence  of  noise,  it  can  be  benef  cial  to  f  Iter 
out  specif  c  ranges  of  wavenumbers  —  ones  that  correspond  to  non-SRBR  modes 
—  because  they  tend  to  cause  striations  whose  slopes  are  unpredictable  and  thus 
corrupt  the  striations  that  are  useful. 

The  f  ltering/beamforming  concepts  in  this  chapter  were  discussed  in  the  con¬ 
text  of  passive  sensing  with  an  array,  but  they  could  likely  be  applied  to  active  arrays 
in  the  context  of  TRM  or  active  sonar,  similar  to  the  ideas  discussed  in  [40].  For 
example,  in  the  context  of  variable  range  TRM,  one  could  have  the  array  only  emit 
sound  over  the  range  of  wavenumbers  corresponding  SRBR  modes  which  typically 
have  Brnj  «  1,  thus  improving  the  ability  to  predict  how  the  range  of  the  focal  spot 
will  change  with  frequency. 
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Chapter  6 

A  relationship  between  the 
waveguide  invariant  and 
wavenumber  integration 


Chapters  3-5  analyzed  the  waveguide  invariant  and  its  applications  using  normal 
modes.  It  was  pointed  out  in  Sec.  2.3  that  the  waveguide  invariant  can  also  be  an¬ 
alyzed  using  ray  theory.  This  chapter  presents  original  research  that  attempts  to 
complete  the  ray/mode/wavenumber-integration  ‘trifecta.’  1  Although  the  research 
presented  in  this  chapter  is  mostly  of  theoretical  interest,  but  may  have  some  prac¬ 
tical  use  in  waveguides  with  attenuating  elastic  media. 


6.1  Introduction 

Although  the  waveguide  invariant  is  usually  derived  and  interpreted  using  nor¬ 
mal  modes  [11],  ray  theory  can  also  be  used  to  interpret  the  waveguide  invariant 
[18,  24,  14],  The  present  analysis  relates  the  waveguide  invariant  to  another  com¬ 
mon  method  for  calculating  the  acoustic  f  eld  in  a  waveguide:  wavenumber  inte¬ 
gration.  It  will  be  shown  that  in  some  cases,  the  waveguide  invariant  can  be  “seen” 

'This  chapter  is  based  on  “A  relationship  between  the  waveguide  invariant  and  wavenumber  inte¬ 
gration”  by  Kevin  L.  Cockrell  and  Henrik  Schmidt  [The  Journal  of  the  Acoustical  Society  of  America 
Express  Letters.  July,  2010]. 
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in  a  plot  of  the  autocorrelation  of  the  wavenumber-integration  kernel.  This  result 
is  somewhat  intuitive  because  the  autocorrelation  of  the  wavenumber-integration 
kernel  represents  horizontal  wavenumber  differences,  which  is  what  the  waveguide 
invariant  is  def  ned  in  terms  of  when  using  a  normal  mode  description  of  the  acous¬ 
tic  f  eld. 

The  present  analysis  is  for  a  range  independent  waveguide  with  planar  geome¬ 
try  so  that  the  wavenumber  transform  is  a  Fourier  transform  (as  opposed  to  a  Bessel 
transform  for  cylindrical  geometry).  Because  the  usual  derivation  of  the  waveguide 
invariant  ignores  cylindrical  spreading,  assuming  planar  geometry  is  not  a  limita¬ 
tion. 


6.2  The  waveguide  invariant  and  normal  modes  in  a  pla¬ 
nar  waveguide 


A  derivation  of  the  waveguide  invariant  using  normal  modes  was  given  in  Sec.  2.3 
(originally  [18]  and  [11,  Sec.  6.7.2]),  using  apoint  source  (cylindrical  geometry)  in 
a  waveguide.  Here  we  present  a  similar  derivation,  but  using  a  planar  geometry.  For 
planar  geometry,  the  complex  pressure  can  be  written  as  as  shown  in  [33,  Eq.  5.26]. 

eikXm(u)\x\ 

P(X,Z,U)  (X  22^m{zs)^m{z)— - -  (6.1) 

m  Kxm\U) 

where  x  is  the  coordinate  pointing  directly  away  from  the  source  and  kxrn(uj )  is  the 
horizontal  wavenumber  for  mode  m  at  a  temporal  frequency  of  to.  Analogous  to 
the  derivation  in  Sec.  2.3,  defne 

Bm  =  T-  (6.2) 

kxm  (W ) 
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The  scalar  acoustic  intensity  is  then  the  pressure  times  its  complex  conjugate: 


I(x,Uj)  =p(x,Uj)  -p(x,Lu) 


oc 


EB1+2E  cos(AA;m/ (cu)x) 

1  m.^l 


(6.3) 

(6.4) 


where  the  overline  indicates  a  complex  conjugate  and  A kmi(uj)  =  kxm(u)  — 
kxi(u).  Eq.  (6.4)  shows  that  the  intensity  at  a  fxed  lo  is  a  sum  of  cosines,  each 
of  which  has  a  spatial  frequency  in  the  x  coordinate  that  depends  on  the  difference 
between  a  pair  of  modes’  horizontal  wavenumbers. 

It  was  shown  in  Sec.  2. 3. 2  that  in  an  ideal  waveguide,  /3  «  1  because  A kmi(oj)  oc 
1/cu  for  modes  far  from  cutoff.  That  implies  that  the  x-coordinate  spatial  frequen¬ 
cies  in  the  intensity  depend  on  u  like  1/oj  [26]. 


6.3  Relating  the  waveguide  invariant  to  wavenumber  in¬ 
tegration 

In  the  previous  section  it  was  shown  that  the  x-coordinate  spatial  frequencies  in 
the  acoustic  intensity  are  determined  by  horizontal  wavenumber  differences.  In 
this  section,  we  relate  the  wavenumber  integration  kernel  to  the  x-coordinate  spa¬ 
tial  frequencies  in  the  intensity  I(x,  u),  which  reveals  a  relationship  between  the 
wavenumber  integration  kernel  and  the  waveguide  invariant. 

To  obtain  the  complex  pressure  p(x,u)  using  wavenumber  integration,  the 
Helmholtz  equation  for  p(x,  uS)  is  transformed  into  the  wavenumber  domain  us¬ 
ing  the  Fourier  transform  pair  [33,  Eqs.  2.85  and  2.86]: 

/OO 

^{kx)eik*xdkx  (6.5) 

-OO 

^(kx)  =  —  /  p(x)e~lkxXdx  (6.6) 

J  —  OO 
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The  boundary  conditions  are  then  written  as  a  function  of  kx,  and  the  resulting 
equations  are  solved  in  the  kx  wavenumber  domain  at  a  single  temporal  frequency 
to  yield  the  wavenumber  kernel  'S?(kx,uj),  which  represents  the  magnitude  and 
phase  of  the  spatial  frequency  components  kx  of  the  complex  pressure  f  eld  at  a  par¬ 
ticular  temporal  frequency  uj.  T  (kx .  uj  )  is  then  transformed  (back)  to  the  x  domain 
using  Eq.  (6.5),  yielding  p(x,u).  The  scalar  acoustic  intensity  can  then  be  com¬ 
puted  by  multiplying  the  pressure  by  its  complex  conjugate:  I ( x ,  uj)  =  p(x,  uj)  ■  p(x,  uj). 
Note  that  the  wavenumber  kernel  T  ( kx ,  uj)  for  a  cylindrical  geometry  is  identical  to 
that  for  a  planar  geometry;  only  the  integral  transform  used  to  obtain  p(x,  uj)  from 
'I >(x,uj )  is  different. 


In  order  to  relate  the  waveguide  invariant  to  the  wavenumber  integration  ker¬ 
nel  4 ’(kx,u ),  we  seek  to  establish  a  relationship  between  '■\)(kx,uj)  and  the  x- 
coordinate  spatial  frequencies  in  the  intensity  I (x,  uj).  This  relationship  is  provided 
by  the  Wiener-Khinchin  Theorem  [63,  64],  which  will  now  be  used  to  show  that  the 
magnitude  of  x-coordinate  spatial  Fourier  Transform  of  I(x,uj)  can  be  calculated 
from  the  autocorrelation  of  the  wavenumber  kernel.  Note  that  although  the  Wiener- 
Khinchin  Theorem  is  well-known  for  its  use  in  statistical  spectrum  estimation,  our 
application  is  purely  deterministic. 


The  autocorrelation  of  the  wavenumber  kernel  is: 

/OO 

^(kx)^(kx  +  A kx)  dkx  (6.7) 

-OO 

where  the  overline  indicates  a  complex  conjugate.  Following  the  standard  deriva¬ 
tion  of  the  Wiener-Khinchin  Theorem,  insert  Eq.  (6.6)  into  Eq.  (6.7)  and  perform  a 
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series  of  algebraic  manipulations  [63,  64]: 


/OO  /  1  poo 

U -J  p(xi)  eikxXl  dx\ 

x  J  p{x 2)  e~l(kx+Ak*)x2  dx^j  dkx  (6.8) 

1  poo  /*oo  /*oo 

=  /  /  /  p(xi)elkxX1p(x2)e~t('kx+Akx'>X2dkxdxidx2 

J  —  OO  J  —  OO  J  —  OO 

1  /*oo  POO  POO 

— 2  /  /  /  p{x\)  p{x2)  elkx^Xl~X2>  e~lAkxX2  dkx  dx\  dx2 

J  — OO  J  —  OO  J  —  OO 


47T2 

1 


'  — OO  «/  — OO  t/  — OO 
poo  POO 

=  ~r?  I  /  p(.Ti)  h(xi  -  x2)p(x2)  e~lAkxX2  dxi  dx2 

J  —  00  «/  —  OO 

1  r00 

=  -^2/  p(xi)p(.Ti)e"*Afca!3:i  dxi 


r  —  OO 
POO 


=  ^LJ{xl)e 


—iAkxxi  , 


(6.9) 


Eq.  (6.9)  shows  that  the  autocorrelation  of  the  wavenumber  kernel  is  propor¬ 
tional  to  the  x-coordinate  spatial  Fourier  transform  of  the  scalar  acoustic  intensity 
/(x).  If  one  were  to  decompose  /(x)  into  all  of  its  spatial  frequencies  compo¬ 
nents,  the  relative  magnitude  of  those  spatial  frequencies  could  be  determined  by 
the  autocorrelation  of  the  wavenumber  kernel. 

This  result  can  be  related  to  the  normal  modes  description  of  the  acoustic  in¬ 
tensity  given  in  Eq.  (6.4).  The  wavenumber  kernel  4 >{kx )  has  peaks  at  values  of 
kx  corresponding  to  the  modal  horizontal  wavenumbers.  Thus  the  autocorrela¬ 
tion  of  T  ( kx )  will  be  large  at  A kx  values  corresponding  to  the  differences  in  the 
modal  horizontal  wavenumbers  —  precisely  the  x-coordinate  spatial  frequencies  of 
/ (x,  uj)  shown  in  Eq.  (6.4). 

The  peaks  (local  maxima)  of  the  autocorrelation  of  the  wavenumber  kernel  cor¬ 
respond  to  the  modal  horizontal  wavenumber  differences  ( Akmi(u )),  so  the  peaks’ 
dependence  on  frequency  will  be  the  same  as  Akmi (cu)’s  dependence  on  frequency. 
Because  the  waveguide  invariant  describes  Ak’m^(cu)’s  dependence  on  frequency, 
the  waveguide  invariant  also  describes  how  the  peaks  of  the  autocorrelation  of  the 
wavenumber  kernel  depend  on  frequency.  The  next  section  analyzes  this  further. 
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6.4  An  ideal  example 


The  waveguide  invariant  is  well  understood  for  an  ideal  waveguide  because  the 
horizontal  wavenumbers  can  be  calculated  analytically.  For  mode  pairs  in  an  ideal 
waveguide  where  both  modes  are  far  from  cutoff,  A kmi(ui)  is  approximately  pro¬ 
portional  to  1/cu,  which  corresponds  to  (5  «  1  (see  Eq.  (2.22)  or  [26]).  So  one 
would  expect  the  A kx  location  of  the  peaks  in  C(Akx,uj )  for  an  ideal  waveguide 
to  depend  on  u  like  l/u.  We  now  show  that  this  is  indeed  the  case. 

The  wavenumber  kernel  for  an  ideal  waveguide  is  [33,  Eq.  2.143] 


^(kx,u)  oc 


sin kzz  sin kz(D— zs) 
kz  sin  kzD 
sin  kzzs  sin  kz(D—z) 
kz  sin  kzD 


Z  <  Zs 
z>  zs  . 


(6.10) 


where  kz  =  \J (^)2  —  kx  and  D  is  the  depth  of  the  waveguide.  if!(kx,u)  has 
poles  when  kzD  =  mir  for  positive  integers  m,  or  equivalently  at  kx  values  cor¬ 
responding  to  the  modal  horizontal  wavenumbers.  Those  poles  will  then  depend 
on  frequency  in  the  same  manner  as  the  modal  horizontal  wavenumbers  will,  and 
so  the  Akx  location  of  the  peaks  of  C(Akx,uj )  will  depend  on  uo  in  the  manner 
predicted  by  the  waveguide  invariant  (oc  1/w). 

To  visualize  this  dependence,  a  plot  of  '£(kx,u)  is  shown  in  Fig.  6.1(a)  along 
with  its  autocorrelation  in  Fig.  6. 1  (b).  The  black  lines  in  Fig.  6. 1  (b)  correspond  to  a 
few  potential  striation  paths  predicted  by  the  waveguide  invariant  with  (3  =  1  (lines 
with  Akx  oc  1/oj).  It  can  be  seen  that  there  are  striations  in  Fig.  6.1(b)  that  do 
not  match  up  well  with  the  black  lines.  This  discrepancy  appears  to  contradict  the 
analysis  in  the  previous  section.  But  if  we  remember  that  /3  is  only  approximately 
equal  to  1  for  mode  pairs  where  both  modes  are  far  from  cutoff,  or  equivalently 
when  kx  is  close  to  k,  then  the  discrepancy  makes  sense  because  the  wavenumber 
kernel  contains  horizontal  wavenumber  components  that  are  close  to  cutoff  and 
thus  are  not  well  described  by  j3  =  1.  To  address  this  issue,  Fig.  6.1(c)  shows  the 
wavenumber  kernel  using  only  kx  values  close  to  k  (specif  cally,  |  k  <  kx  <  k)  and 
Fig.  6.1(d)  shows  the  resulting  autocorrelation  along  with  black  lines  corresponding 
to  [3  =  1.  The  striations  in  Fig.  6.1(d)  match  very  well  with  the  black  lines  because 
only  horizontal  wavenumbers  that  are  far  from  cutoff  are  included. 
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Figure  6.1:  All  subfgures  are  for  a  100  m  deep  ideal  waveguide  with  zs  =  16  m 
and  z  =  73m,  and  are  in  units  of  dB  re:  an  arbitrary  reference,  (a)  \'-\l(kx-co)\, 
the  wavenumber  kernel  as  a  function  of  horizontal  wavenumber  and  frequency,  (b) 
\C(Akx,uj)\,  the  autocorrelation  of  wavenumber  kernel  shown  in  subfgure  (a)  as 
a  function  of  wavenumber  difference  A kx  and  frequency  oj.  The  black  lines  are 
example  striation  paths  corresponding  3  =  1  (lines  with  A kx  oc  1  /oj).  (c)  Same 
as  subfgure  (a),  but  only  including  kx  values  far  from  cutoff  (| k  <  kx  <  k).  (d) 
\C(Akx,  Co-)  | ,  the  autocorrelation  of  the  wavenumber  kernel  shown  in  subf  gure  (c). 
The  black  lines  are  example  striation  paths  corresponding  3  =  1  (lines  with  A kx  oc 
1/c o),  and  match  the  actual  striation  paths  more  closely  than  those  in  subf  gure  (b) 
because  subf  gure  (d)  only  includes  kx  values  far  from  cutoff. 
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6.5  Relevance 


Many  applications  of  the  waveguide  invariant,  such  as  the  one  in  Ch.  3,  require  an 
assumption  about  the  value  of  f3.  It  is  often  correct  to  assume  that  /3  «  1  in  shallow- 
water  waveguides,  but  this  assumption  is  not  always  valid,  as  discussed  earlier 
in  this  thesis  and  in  [18].  Consequently,  numerical  modeling  is  sometimes  used 
to  determine  the  approximate  value  of  (3  in  a  given  environment  [51],  Although 
this  can  be  done  by  simulating  the  acoustic  f  eld  itself,  more  insight  can  often  be 
gained  by  calculating  more  fundamental  quantities  such  as  the  modal  horizontal 
wavenumbers  or  the  ray  parameters  (horizontal  slowness,  cycle  distance,  etc.  -  see 
[24,  Eqs.  11-15]),  which  can  be  related  to  the  value  of  the  waveguide  invariant. 

However,  normal  modes  and  ray  theory  are  not  conducive  to  describing  the 
acoustic  f  eld  in  some  environments,  such  as  those  involving  attenuating  elastic 
media  (especially  if  one  is  interested  in  the  acoustic  f  eld  inside  of  the  elastic  me¬ 
dia).  In  those  cases,  understanding  how  the  wavenumber-integration  kernel  relates 
to  the  waveguide  invariant  may  allow  one  to  gain  insights  that  would  be  diff  cult  to 
obtain  otherwise. 

For  example,  the  concepts  described  in  this  chapter  may  be  useful  for  study¬ 
ing  the  effect  that  surface  ice  has  on  the  value  of  j3.  The  striations  in  a  plot  like 
Fig.  6.1(d),  but  generated  for  an  ice-covered  waveguide  rather  than  an  ideal  one, 
could  reveal  whether  it’s  reasonable  to  assume  j3  =  1  in  such  an  environment. 
More  specif  cally,  one  could  use  a  plot  like  Fig.  6.1(d)  to  determine  if  hydrophones 
sitting  on  surface  ice  (or  embedded  in  the  seaf  oor)  would  record  the  the  same 
range-frequency  waveguide  invariant  striations  that  a  hydrophone  in  the  water  col¬ 
umn  would  record. 

6.6  Chapter  summary 

The  waveguide  invariant  implies  a  specif  c  dependence  of  the  acoustic  intensity’s 
x-coordinate  spatial  frequencies  on  the  temporal  frequency  ui.  When  the  geom¬ 
etry  of  the  problem  is  planar,  the  wavenumber  integration  transform  is  a  Fourier 
Transform,  and  so  the  Wiener-Khinchin  Theorem  can  be  used  to  relate  the  auto¬ 
correlation  of  the  wavenumber  kernel  to  the  x-coordinate  spatial  frequencies  in  the 
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acoustic  intensity.  A  2-d  plot  of  the  autocorrelation  of  the  wavenumber  kernels  ver¬ 
sus  “wavenumber  lag”  A kx  and  temporal  frequency  u>  exhibits  striation  patterns 
that  can  be  explained  by  the  waveguide  invariant. 

The  relationship  between  the  waveguide  invariant  and  wavenumber  integration 
is  not  as  direct  as  the  waveguide  invariant’s  relationship  to  normal  modes  or  ray 
theory.  However,  the  relationship  described  in  the  present  analysis  may  be  useful 
in  certain  situations  because  wavenumber  integration  can  be  used  in  waveguides 
that  are  diff  cult  to  model  using  normal  modes  or  ray  theory. 
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Chapter  7 

Conclusion 


The  waveguide  invariant  8  is  both  embarrassingly  simple  (why  did  it  take  so  long 
to  explain  the  striations  in  spectrograms!?),  and  fantastically  complicated  (hidden 
behind  so  many  layers  of  mathematical  manipulation,  it’s  a  wonder  it  was  ever 
discovered!). 

To  explain  waveguide  invariant  range-frequency  striations,  one  must  transform 
the  time-domain  wave  equation  into  the  frequency  domain.  Then  transform  that 
frequency-domain  equation  into  the  horizontal  wavenumber  domain,  and  then  sepa¬ 
rate  it  into  a  depth-dependent  equation  and  a  range-dependent  equation.  The  range- 
dependent  equation  provides  the  elkrmr  terms  that  interfere  with  one  another.  The 
depth-dependent  equation,  which  often  cannot  be  solved  analytically,  provides  the 
krm  which  depend  on  u>.  The  difference  between  two  krm s,  and  how  that  differ¬ 
ence  depends  on  c o,  determines  the  waveguide  invariant  f3.  And  8  determines  the 
slopes  of  the  striations  in  a  spectrogram  from  a  hydrophone  being  towed  away  from 
a  source. 

Small  changes  in  the  SSP  and  in  the  seaf  oor  properties  can  cause  large  changes 
in  the  structure  of  the  acoustic  f  eld.  So  it  is  somewhat  surprising  that  a  single  scalar 
value  of  8  accurately  describes  the  range-frequency  striations  that  appears  in  the 
acoustic  f  eld  under  such  a  wide  range  of  SSPs  and  seaf  oors. 

This  thesis  built  upon  previous  research  to  understand  why  8  is  such  a  robust 
description  of  the  acoustic  f  eld,  and  to  develop  fundamental  signal  processing  tech¬ 
niques  for  the  waveguide  invariant. 
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This  chapter  summarizes  this  thesis  and  discusses  its  signif  cance.  Although 
this  thesis  made  signif  cant  contributions,  much  is  left  to  be  understood;  several 
suggestions  for  future  research  are  given  in  this  chapter. 


7.1  Thesis  summary 


Chapter  3  assumed  j3  =  1,  which  allowed  for  passive  range  estimation.  Although 
similar  research  had  already  been  performed,  Chapter  3  quantif  ed  the  signal  pro¬ 
cessing  techniques  and  derived  some  fundamental  quantities  such  as  the  minimum 
source  bandwidth  required  to  observe  range-frequency  striations.  Chapter  3  also 
provided  an  experimental  example  showing  that  /3  «  1  in  a  realistic  waveguide 
with  a  non-uniform  SSP  and  a  seaf  oor  with  unknown  acoustic  properties. 

Chapter  4  analyzed  why  (3  «  1  in  shallow-water  waveguides,  and  how  the 
SSP  and  seaf  oor  affect  j3.  The  original  paper  on  the  waveguide  invariant  gave 
some  brief  hand-waiving,  albeit  substantive,  arguments  for  this,  and  solved  for  /3 
analytically  in  a  few  special  cases.  Other  researchers  investigated  f3  numerically, 
and  analytically  using  ray  theory.  In  Ch.  4,  a  method  was  derived  using  normal 
modes  that  provided  not  only  a  quantitative  relationship  between  the  SSP  and  (3, 
but  also  a  relationship  that  could  also  be  understood  intuitively  through  the  function 
c t>i{Sp ).  This  method  was  used  to  prove  that  f3  >  0  for  all  surface-ref  ected  bottom- 
ref  ected  (SRBR)  modes  in  a  vacuum-bounded  waveguide  regardless  of  the  SSP.  It 
was  also  used  for  a  qualitative  argument  suggesting  that  13  ~  1  for  SRBR  modes  in 
a  vacuum-bounded  waveguide,  regardless  of  the  SSP.  (Unfortunately  the  analysis 
in  chapter  4  did  not  prove  that  ~  1  for  SRBR  modes.) 

Chapter  5  drew  from  the  understanding  gained  in  Ch.  4,  and  showed  how  an 
array  can  be  used  to  not  only  f  Iter  out  noise,  but  to  flter  out  non-SRBR  modes 
leaving  only  the  SRBR  modes  which  are  likely  have  (3mi  ~  1. 

Chapter  6  related  the  waveguide  invariant  to  wavenumber  integration,  develop¬ 
ing  another  tool  to  study  the  waveguide  invariant. 
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7.2  Thesis  Significance 


Applications :  The  waveguide  invariant  has  been  applied  to  a  broad  range  of  prob¬ 
lems  in  ocean  acoustics:  passive  sonar,  active  sonar,  array  processing,  time-reversal 
mirrors,  and  more.  (See  Appendix  A  for  a  review  of  the  waveguide  invariant  lit¬ 
erature).  Many  of  those  applications  could  benef  t  from  the  concepts  discussed  in 
this  thesis.  For  example,  a  transmitting  TRM  array  could  emit  sound  primarily  at 
horizontal  wavenumbers  that  correspond  to  SRBR  modes,  ensuring  that  3  ~  1  and 
thus  increasing  the  predictability  of  the  focal  spot  shift  with  frequency  (details  in 
next  section). 

Fundamental  Understanding :  The  research  in  this  thesis  was  applied  exclu¬ 
sively  to  range  independent  waveguides.  Flowever,  the  waveguide  invariant  is  com¬ 
monly  applied  to  range  dependent  and  azimuthally  dependent  waveguides.  The 
method  for  calculating  3  in  Ch.  4  could  be  applied  to  understanding  how  range- 
frequency  striations  behave  in  such  environments.  It  could  also  be  used  for  under¬ 
standing  how  the  range-frequency  striations  change  in  time  due  to  oceanographic 
events  like  internal  waves  (details  in  next  section). 

Also,  the  analysis  in  Ch.  6  relating  the  waveguide  invariant  to  wavenumber 
integration  may  be  useful  for  studying  the  waveguide  invariant  under  conditions 
which  are  not  conducive  to  ray  theory  or  normal  modes,  such  as  the  acoustic  f  eld 
in  a  lossy  elastic  seaf  oor. 

7.3  Specific  suggestions  for  future  research 

This  section  is  a  list  of  suggestions  for  future  research.  Some  suggestions  are  nat¬ 
ural  continuations  of  the  research  presented  in  this  thesis,  and  some  are  simply 
“holes”  in  the  current  understanding  of  the  waveguide  invariant.  The  diff  culty  of 
executing  these  suggestions  ranges  from  “assign  as  a  homework  problem”  to  “write 
a  multi-year  multi-institution  research  grant  proposal.” 

Range  dependent  environments 

The  waveguide  invariant  3  can  be  def  ned  for  mildly  range-dependent  environments 
using  the  adiabatic  approximation  [18].  D’Spain  et  al.  extended  that  to  range-  and 
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azimuthally-dependent  environments  using  the  “N  by  2D”  approximation  [21],  The 
method  for  calculating  (3  presented  in  Ch.  4  could  be  applied  to  (or  even  simply 
inserted  into)  [21,  Eq.  14  and/or  18]. 

The  concept  of  the  waveguide  invariant  may  even  be  able  to  be  extended  to 
strongly  range-dependent  environments  where  mode  coupling  occurs.  [53]  has  an 
expression  for  how  the  modal  amplitudes  Arn  change  as  a  function  of  range  due  to 
mode-coupling.  One  may  be  able  to  insert  that  expression,  using  the  mode  shapes 
as  calculated  from  the  WKB  approximation,  into  the  derivation  of  the  waveguide 
invariant  shown  in  Sec.  2.3  to  see  how  much  mode  couping  can  occur  before  the 
range-frequency  striations  are  destroyed. 

Forcing  /3  «  1  without  an  array  (AUV  environmental  adaptation) 

Chapter  4  showed  that  surface-ref  ected  bottom-ref  ected  (SRBR)  modes  tend  to 
have  f3  «  1.  Chapter  5  showed  that  if  one  wants  to  observe  only  SRBR  modes, 
one  can  use  an  array  to  f  Iter  out  certain  wavenumbers.  An  alternative  method  for 
observing  the  striations  from  only  SRBR  modes  (which  tend  to  have  (3  ~  1)  is 
to  place  the  acoustic  receiver  at  the  depth  of  the  sound  speed  maximum.  Doing 
so  will  minimize  the  contribution  to  I(r,  u)  from  non-SRBR  modes,  because  non- 
SRBR  modes  will  have  small  'ip(z)  at  that  depth.  This  could  be  useful  for  any 
application  where  one  would  like  to  assume  (3  «  1,  like  passive  range  estimation. 
An  Autonomous  Underwater  Vehicle  (AUV)  could  adaptively  determine  the  SSP 
maximum,  and  move  to  that  depth.  An  AUV  with  an  array  could  detect  a  source, 
adaptively  head  directly  towards  it,  and  move  to  depth  of  the  SSP  maximum. 

Tighter  bounds  on  j3  for  SRBR  modes 

In  Sec.  4.4,  it  was  mathematically  proven  that  fi(Sp)  >  0  for  SRBR  modes  in  a 
vacuum-bounded  waveguide,  regardless  of  the  SSP.  The  examples  in  that  chapter 
also  showed  that  typically  (3(SP)  ~  1  for  SRBR  modes.  A  “hand-waiving”  argu¬ 
ment  was  used  to  explain  this:  SRBR  modes  behave  similarly  to  modes  in  an  ideal 
waveguide,  which  have  /3  «  1  (see  last  paragraph  of  Sec.  4.2).  By  using  the  equa¬ 
tions  in  Sec.  4.4,  it  may  be  possible  derive  tighter  bounds  on  (3,  or  to  quantitatively 
explain  why  (3  «  1  for  SRBR  modes  regardless  of  the  SSP. 
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Waveguide  invariant  array  processing  for  active  transmissions 

The  array  processing  concepts  presented  in  Ch.  5,  where  one  f  ltered  out  non-SRBR 
which  typically  do  not  correspond  to  (3  «  1,  could  be  applied  to  many  of  the  waveg¬ 
uide  invariant  applications  discussed  in  Appendix  A.  Specif  cally,  transmitting  ar¬ 
rays  in  TRM  may  be  able  to  apply  similar  concepts,  increasing  the  ability  accurate 
predict  how  the  range  of  the  focal  spot  will  shift  with  frequency.  The  work  pre¬ 
sented  in  [49]  may  also  benef  t  from  such  f  ltering,  although  the  mode  coupling 
may  subdue  the  desired  effect. 

Waveguide  invariant  array  processing  in  non-uniformly  spaced  arrays 

The  array  processing  (wavenumber  f  ltering)  method  in  Ch.  5  was  based  on  a  stan¬ 
dard  time-domain  f  Iter  design  technique  (f  i  r  1  s  in  Matlab).  Because  time  domain 
f  Iter  design  techniques  usually  assume  the  samples  are  spaced  uniformly  in  time, 
they  cannot  be  applied  to  arrays  where  the  elements  are  non-uniformly  spaced. 
If  one  wanted  to  do  wavenumber  f  ltering  similar  to  that  in  Ch.  5,  but  for  a  non- 
uniformly  spaced  array,  one  would  have  to  use  (or  develop)  f  Iter  design  techniques 
that  can  handle  non-uniformly  spaced  arrays.  The  methods  in  [61,  Ch.  3]  may  be 
useful  for  this. 


l/r  spreading 

Quantify  the  effect  of  the  l/r  spreading  that  is  ignored  when  deriving  the  waveg¬ 
uide  invariant,  in  terms  of  how  it  affects  the  “observed44  value  of  /j  in  the  striation 
slopes  in  I (r,  u).  Note  that  one  can  factor  out  the  1  /  yfr  in  the  Am  and  .4/  terms  in 
Eq.  2.10. 

Including  attenuation  when  deriving  f3 

Re-derive  the  waveguide  invariant  but  allow  for  the  horizontal  wavenumbers  to  be 
complex  so  that  they  can  take  into  account  bottom  loss.  The  result  will  likely  be 
like  Eq.  2.10,  but  with  an  exponential  decay  tacked  on  to  each  cosine  term.  The 
footnote  in  Sec.  2.3.1  may  help  with  this. 
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Group  speeds  and  phase  speeds 


[10,  Eq  4.5.15]  gives  an  integral  expression  for  the  group  speed  times  the  phase 
speed,  um  ■  vm,  for  a  c(z)  that  is  constant  above  some  depth  but  can  vary  with  z 
below  that  depth.  That  expression  may  be  able  to  be  used  to  calculate  /3,  without 
even  using  the  WKB  approximation.  Or  it  may  provide  some  general  bounds  or 
insights  on  the  waveguide  invariant  (in  that  specif  c  environment). 


Striation  slopes  change  with  range  and  frequency 

As  discussed  in  Sec.  2.2.3  and  Sec.  2.3.3,  the  waveguide  invariant  (3mi  is  typi¬ 
cally  def  ned  using  adjacent  modes  (the  Taylor  series  expansion  of  group  and  phase 
slowness  plotter  versus  mode  number).  The  justif  cation  for  this  is  based  on  a  prop¬ 
erty  of  the  modal  sum  discussed  in  Sec.  2.2.3.  The  group  (or  groups)  of  modes 
that  dominate  the  modal  sum  is  determined  by  the  range  and  frequency  of  the  re¬ 
ceiver.  Thus,  the  range-frequency  striation  slopes  (locally  in  the  (r,  /)  plane)  will 
be  different  at  different  ranges  and  frequencies,  as  is  observed  in  Fig.  2.6(c).  This 
phenomena  is  not  observed  as  strongly  in  shallow-water  waveguides,  presumably 
because  shallow-water  waveguides  have  many  fewer  modes  and  less  of  those  modes 
are  non-SRBR  modes. 

It  would  be  worthwhile  to  quantitatively  investigate  when  the  argument  in  Sec.  2.2.3 
applies,  and  how  it  could  be  used  to  predict  which  group  of  modes  (and  thus  which 
value  of  /3)  is  relevant  at  a  given  range  frequency.  Perhaps  there  are  certain  ranges 
and  frequencies  where  one  can  assume  (3  «  1  in  a  deep-water  waveguides,  as  sug¬ 
gested  by  some  of  the  striations  around  60  km  in  Fig.  2.6(c)  (which  happens  to 
be  the  convergence  zone  distance...).  This  investigation  could  be  done  using  the 
argument  in  Sec.  2.2.3  —  looking  at  kr  as  a  function  of  m  and  u.  In  addition,  ray 
theory  may  also  provide  good  insights  into  this,  especially  as  it  is  portrayed  in  [20, 

Fig.  3.8]  (which  is  the  original  reference  for  Sec.  2.2.3). 

Another  approach  for  studying  this  would  be  to  investigate  the  relationship  be¬ 
tween  the  ( AmAi  cos(A kmir))  terms  in  Eq.  2.10,  and  the  argument  in  Sec.  2.2.3. 
Figure  2.2  shows  the  contribution  to  the  complex  pressure  from  each  mode,  but 
each  ( AmAi  cos(A  kmir))  term  represents  the  contribution  to  the  magnitude  of  the 
pressure  from  2  modes. 
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Upper  limits  on  range  and  frequency  for  observing  striations 

Waveguide  invariant  striations  occur  when  the  environment  is  not  strongly  range- 
and/or  azimuthally-dependent  (i.e.,  no  mode  coupling),  and  the  source  is  broad¬ 
band.  As  the  frequency  and  range  increase,  scattering  from  inhomogeneities  in 
the  environment  (boundary  roughness,  SSP  variations)  tends  to  smear  out  or  com¬ 
pletely  destroy  the  waveguide  invariant  striations.  This  effect  can  be  seen  in  Fig.  2.5, 
where  the  striation  widths  get  wider  as  the  range  increases  (contrary  to  what  Eqs.  (3.6) 
and  (3.5)  predict  for  a  perfectly  homogeneous  waveguide).  Rouseff  and  Kuz’kin 
have  done  some  research  on  this  topic  [52,  42,  41],  but  there  is  certainly  more  re¬ 
search  that  could  be  done  with  relative  ease  given  the  state  of  understanding  of 
acoustic  scattering  in  the  ocean  waveguide. 

The  upper  limit  in  range  and  frequency  at  which  the  striations  would  appear 
will  depend  on  the  inhomogeneities  of  a  particular  environment.  Such  upper  limits 
could  obviously  be  measured  experimentally  by  towing  a  hydrophone  away  from 
a  very  broadband  source.  In  addition,  the  literature  in  Appendix  A  provides  much 
experimental  data. 

On  the  theoretical  side,  there  is  a  large  body  of  literature  studying  scattering 
in  the  ocean;  surely  it  could  be  applied  to  determining  how  far  away  and  how 
high  frequency  waveguide  invariant  striations  could  be  observed.  To  frst  order, 
surface  and  seaf  oor  roughness  act  simply  as  a  loss-mechanism  [33,  Sec.  1.7].  More 
sophisticated  models  would  likely  need  to  be  used.  [60,  Ch.  6]  may  provide  a  good 
starting  point,  as  would  [36,  Chs.  4  and  6]. 

Another  approach  might  be  to  use  the  expressions  used  to  generate  [16,  Fig.  5(a)], 
which  could  be  analyzed  as  a  function  range  and  frequency.  [2,  17]  may  also  pro¬ 
vide  useful  insights. 
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Appendix  A 

Review  of  waveguide  invariant 
literature 


This  appendix  is  a  list  of  references  about  the  waveguide  invariant,  with  each  ref¬ 
erences  loosely  categorized  based  on  its  application.  Nearly  all  articles  in  The 
Journal  of  the  Acoustical  Society  of  America  and  IEEE  Oceanic  Engineering  relat¬ 
ing  to  the  waveguide  invariant  are  included  here;  some  of  the  Russian  literature  is 
also  included. 1  [47]  gives  a  brief  summary  of  the  Russian  literature  on  the  waveg¬ 
uide  invariant  and  similar  phenomena,  including  several  of  references  which  are 
not  reviewed  in  this  appendix. 

A.l  General  review 

Fundamental  research  on  the  waveguide  invariant 

The  f  rst  paper  on  the  waveguide  invariant  is  [1 8]  by  Chuprov.  This  paper  discusses 
the  waveguide  invariant  from  the  perspective  of  normal  modes  and  ray  theory.  It 
focuses  primarily  on  range-frequency  striations,  but  also  discusses  some  changes 
with  receiver  depth.  It  focuses  mostly  on  deep-water  propagation  but  discusses 
shallow-water  waveguides  as  well.  It  analytically  calculates  the  waveguide  invari- 

*The  journal  Soviet  Physical  Acoustics  turned  into  Acoustical  Physics  sometime  in  the  1990s. 
Articles  from  Acoustical  Physics  are  available  online  from  Springer;  many  articles  in  Soviet  Physical 
Acoustics  have  been  translated  into  English  and  were  able  to  be  tracked  down  MIT’s  superb  librarians. 
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ant  for  ideal  waveguides  and  surface-trapped  modes  in  an  n2- linear  waveguide.  It 
presents  numerical  and  experimental  data,  and  shows  that  in  most  cases  the  Doppler 
effect  due  to  source/receiver  motion  has  little  effect  on  range-frequency  striations. 
It  also  brief  y  mentions  how  to  account  for  mild  range  dependence  using  the  adia¬ 
batic  approximation  (although  other  references  give  a  more  detailed  derivation).  In 
general,  this  paper  is  dense  and  fairly  diff  cult  to  follow  unless  the  reader  has  a  very 
deep  understanding  of  ocean  acoustics.  It’s  a  ‘must  read’,  but  probably  only  after 
reading  other  simpler  derivations  of  the  waveguide  invariant. 

In  [11,  Sec.  6.7.2],  Brekhovskikh  and  Lysanov  give  a  derivation  and  discus¬ 
sion  of  the  range-frequency  waveguide  invariant  in  somewhat  simpler  terms  than 
Chuprov’s  original  paper  does.  2 

[21]  by  D’  Spain  et  al.  is  a  highly  cited  paper  which  presents  more  thorough 
derivations  of  many  of  the  results  in  Chuprov’s  original  paper.  This  paper  extends 
Chuprov’s  range-dependent  expression  for  6  to  environments  that  change  with  az¬ 
imuth  as  well,  using  the  “N  by  2D”  approximation.  Striations  in  experimental  data 
with  a  frequencies  of  70-170  Hz  are  compared  with  predicted  striation  slopes.  Also 
derived  is  an  expression  for  the  striation  trajectories  in  the  spectrogram  of  a  f  xed 
hydrophone,  when  the  a  source’s  position  is  a  function  of  time.  This  paper  provides 
an  approximate  expression  for  f3  in  an  waveguide  with  a  constant  sound  speed  and 
vacuum  bottom,  but  with  range-dependent  bathymetry.  Also  discussed  are  “gener¬ 
alized”  waveguide  invariants,  which  explain  striations  which  results  from  changes 
in  waveguide  depth. 

The  2nd  edition  of  the  textbook  “Computational  Ocean  Acoustics44,  which  is 
still  in  preparation,  will  contain  a  few  sections  that  discuss  the  waveguide  invariant. 
[34,  Sec.  2.4.6]  gives  a  very  nice  introduction  to  the  waveguide  invariant.  There  is 
a  section  in  [34,  Ch.  3]  that  discusses  how  the  waveguide  invariant  can  be  calcu¬ 
lated  from  acoustic  ray  parameters,  which  is  based  on  the  work  in  [22]  and  [24]. 
Expressions  are  derived  for  calculating  the  waveguide  invariant  without  needing  to 
run  full  ray-tracing  models,  similar  to  the  result  in  Sec.  4  of  this  thesis  which  does 
not  require  solving  for  the  horizontal  wavenumbers.  A  section  in  [34,  Ch.  5]  also 
discusses  the  waveguide  invariant,  and  gives  some  numerical  examples. 

In  [26],  Grachev  shows  that  if  a  single  scalar  value  of  (3  is  an  accurate  descrip- 

2This  author  recalls  noticing  an  typographical/mathematical  error  in  their  derivation,  however. 
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tion  of  the  horizontal  wavenumbers,  then  those  horizontal  wavenumbers  must  have 
a  power-law  dependence  on  frequency.  Grachev  also  discusses  an  analogous  result 
or  a  generalized  waveguide  invariant  that  deals  with  changes  in  waveguide  depth. 
Also  discussed  are  (3  for  surface-trapped  modes  modes  in  n2 -linear  waveguides, 
and  (3  for  waterborne  modes  in  asymmetric  and  symmetric  n2- linear  waveguides. 

In  [42],  Kuz’kin  et  al.  study  range-frequency  striations  in  a  range-dependent 
environment  using  experimental  data.  Perturbation  theory  is  used  to  determine 
how  changes  in  the  SSP  will  affect  changes  in  the  observed  striation  pattern,  and 
how  one  might  be  able  to  use  that  for  acoustic  inversions. 

In  [41],  Kuz’kin  uses  a  statistical  correlation  framework  to  study  range-frequency 
striations  (or  frequency  shifts)  due  to  range-dependent  SSP  perturbations. 

In  [15],  Bumekov  discusses  the  boundaries  (range)  of  validity  of  the  adiabatic 
range-dependent  expression  for  (3.  He  analyzes  (3  in  the  case  of  a  waveguide  with 
a  constant  bottom  slope,  and  he  shows  that  (3  as  a  function  of  range  can  change 
quickly  with  range  even  if  the  environment  changes  slowly  with  range  because 
modes  can  go  from  being  SRBR  to  non-SRBR  as  the  bathymetry  changes. 

In  [44],  Lobanov  and  Petukhov  discuss  that  in  range-dependent  environments, 

(3  can  have  a  “singularity”  (/3  ( r\ )  =  0)  which  implies  the  range-frequency  striations 
will  be  completely  parallel  to  the  frequency  axis  in  the  range-frequency  plane.  They 
point  out  that  in  the  time  domain,  j3{r\)  =  0  means  that  pulse  signals  are  not 
subject  to  dispersion  broadening  with  time  at  the  range  r\.  [44]  investigates  the 
range-frequency  striations  when  only  a  few  modes  are  propagating,  and  notes  that 
even  if  many  modes  are  propagating,  the  higher  order  modes  often  attenuate  away 
with  range  leaving  only  a  few  propagating  modes.  The  effects  of  seismic  (elastic) 
waves  on  the  range-frequency-striations  are  also  brief  y  discussed.  Finally,  there  is 
a  brief  mention  of  how  (3  can  change  with  range  and  frequency. 

Note  that  much  of  the  Russian  literature,  such  as  [44],  makes  statements  that 
only  apply  to  waterborne  modes  without  explicitly  stating  that  the  statements  only 
apply  to  waterborne  modes.  For  example:  “under  the  WKB  approximation,  a  group 
speed  has  a  functional  dependence  on  the  phase  speed  [independent  of  frequency].” 
That  statements  is  only  true  for  waterborne  modes,  and  is  not  true  for  bottom  inter¬ 
acting  modes  (Ch.  4  illustrates  why). 

In  [52],  Rouseff  uses  numerical  simulations  to  study  the  effect  of  shallow-water 
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internal  waves  on  /3.  He  discusses  that  the  value  of  3  depends  on  the  type  of  mode, 
and  that  the  “measured44  value  of  3  (i.e.,  the  striations  slope  seen  in  I (r,  cu))  depends 
on  receiver  depth.  He  numerically  studies  the  effect  of  bottom  attenuation  and 
range  on  the  range-frequency  slopes. 

In  [4],  Baggeroer  discusses  methods  for  plotting  and  understanding  the  waveg¬ 
uide  invariant  in  several  canonical  waveguides,  using  normal  mode  simulations. 
Baggeroer  also  shows  how  the  2D-DFT  can  be  used  to  extract  the  value  or  distri¬ 
bution  of  3  from  I (?’,  w). 

In  [51],  Rouseff  and  Spindel  discuss  how  to  use  a  2D-DFT  to  measure  the  dis¬ 
tribution  of  3  from  7(r,  cu),  and  show  that  the  striations  slopes  sometimes  depends 
on  the  source  and  receiver  depths. 

In  [14],  Brown  et  al.  shows  that  when  the  modal  horizontal  wavenumbers  are 
calculated  from  from  the  WKB  approximation  and  the  limit  of  high  frequency 
is  taken,  3  is  equivalent  to  the  “ray  stability  parameter”  a.  Note  that  there  is 
some  overlap  between  the  work  in  [14]  and  the  work  Ch.  4.  [14]  discusses  range- 
dependent  effects,  including  mode-coupling.  [6]  extends  the  work  of  [14]  and  dis¬ 
cusses  how  3  related  to  “acoustic  beam  dynamics.” 

In  [1],  An  et  al.  discuss  how  a  2D-DFT  can  be  used  to  calculate  the  slope  of  the 
waveguide  invariant  striations.  They  show  a  simulated  plot  of  3  versus  frequency 
and  modal  phase  speed  in  a  Pekeris  waveguide.  They  also  show  how  the  value  of  3 
observed  in  7  (?’,  uj)  depends  on  the  source  and  receiver  depths,  and  they  show  some 
experimental  data  from  the  South  China  Sea  with  frequencies  around  370  Hz  and 
ranges  around  13  km. 

In  [13],  Brooks  et  al.  discusses  a  few  methods  of  extracting  the  striation  slopes 
from  7(r,  w). 

[47],  which  is  a  summary  of  the  Russian  literature  on  the  waveguide  invari¬ 
ant,  is  the  f  rst  chapter  of  a  book  which  contains  several  articles  on  the  waveguide 
invariant,  some  of  which  are  not  described  in  this  appendix. 

Active  sonar 

In  [49],  Quijano  et  al.  discuss  how  the  waveguide  invariant  concept  can  be  applied 
to  bi-static  active  sonar.  An  expression  is  given  for  the  pressure  as  a  function  of 
source/receiver/target  geometry,  including  the  mode  coupling  which  depends  on  the 
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model  of  the  target.  That  expression  illustrates  how  /3  from  passive  measurements 
can  be  related  to  active  sonar.  Some  experimental  data  shows  a  measured  /3active 
of  1.4,  which  is  then  compared  to  the  (passively)  measured  values  of  /3  from  other 
papers. 

In  [27],  Chensong  et  al.  applies  the  concept  of  the  active  sonar  waveguide 
invariant  from  [49]  to  an  extended  Kalman  f  Iter  tracking  algorithm.  The  waveguide 
invariant  is  used  to  restrict  the  possible  state  transitions  of  the  state  vector  in  the 
Kalman  f  Iter,  leading  to  better  tracking  performance. 

In  [25],  Goldhahn  et  al.  uses  the  waveguide  invariant  to  describe  reverberation 
in  active  sonar  returns,  allowing  for  an  improved  ability  for  target  detection. 

Passive  sonar  and  non-adaptive  array  processing 

In  [56],  Tao  et  al.  demonstrate  a  method  for  estimating  the  speed  and  closest-point- 
of-approach  of  a  moving  source.  Tao  shows  that  if  (3  =  1  and  the  source  maintains 
a  constant  speed  and  bearing,  parabolic  striations  will  appears  in  a  spectrogram 
from  single  f  xed  hydrophone.  The  properties  of  the  striations  are  extracted  using 
a  Hough  transform,  and  a  range  and  speed  estimate  are  produced.  The  method  is 
tested  on  experimental  data  from  a  pontoon  boat  in  a  9  m  deep  lake,  with  frequen¬ 
cies  of  about  3-5  kHz  and  ranges  of  a  few  hundred  meters. 

In  [23],  D’Spain  et  al.  discuss  striations  in  a  plot  the  beamformed  acoustic  in¬ 
tensity  from  a  vertical  array.  Striations  in  range,  frequency  and  vertical  wavenum¬ 
ber  are  discussed  assuming  the  vertical  array  is  an  ideal  wavenumber-f  Iter  (as  op¬ 
posed  to  [66]  and  Ch.  3,  which  explicitly  take  into  account  the  array’s  non-ideal 
response).  The  results  are  discussed  in  the  context  of  4  different  SSPs  which  can 
be  treated  analytically.  Methods  for  depth  and  range  estimation  are  discussed  and 
demonstrated  using  experimental  data  from  SWellEx-3  (frequencies  of  300-700  Hz 
and  ranges  of  1.5  to  6  km.) 

In  [66],  Yang  shows  that  the  waveguide  invariant  range-frequency  striations 
appear  in  the  output  of  a  horizontal  line  array  beamformer,  with  some  array  gain. 
He  discusses  some  applications  of  this  and  demonstrates  them  with  simulated  data. 

In  [62],  Turgut  et  al.  estimate  the  ratio  of  a  broadband  acoustic  source’s  speed 
to  its  closest-point-of-approach.  Striation  trajectories  are  derived  for  sources  with 
constants  speeds  and  bearings,  similar  to  those  in  [56],  The  striations  observed  in 
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the  beamformer’s  output  from  experimental  data  (frequencies  of  50-150  Hz)  are 
analyzed  using  Hough  transforms.  The  waveguide  invariant  formulation  from  [21] 
is  used  to  take  into  account  a  (3  that  varies  with  range  and  azimuth.  [62]  also 
investigates  the  possibility  of  jointly  estimating  the  value  of  /j  along  with  the  other 
parameters. 

Chapters  3  and  5  of  this  thesis,  and  their  associated  papers  in  JASA  would  be 
in  this  section. 

In  [59],  Thode  et  al.  localize  an  acoustic  source  by  utilizing  the  waveguide 
invariant  in  a  matched  f  eld  processing  (MFP)  scheme.  Thode  shows  that  the  side- 
lobes  of  the  MFP,  when  plotted  versus  frequency  for  a  broadband  source,  follow 
a  trajectory  described  by  3  which  can  be  used  to  accurately  estimate  the  range  to 
the  acoustic  source  assuming  (3  =  1.  Thode  demonstrates  this  experimentally  for 
shallow  water  with  a  Blue  whale  vocalization.  Thode  also  argues  through  analysis 
and  simulation  that  even  if  (3  is  a  negative  number,  as  it  would  be  in  environments 
like  the  Arctic,  one  can  use  the  sidelobe  behavior  to  estimate  the  source  range.  This 
result  is  interesting  because  most  other  applications  of  the  waveguide  invariant  as¬ 
sume  (3=1. 

In  [58],  Thode  presents  a  method  for  “estimating  the  range  of  an  unknown 
broadband  acoustic  source  in  a  [range-independent]  waveguide,  using  a  vertical 
array  and  a  signal  sample  from  another  broadband  source  a  known  location  relative 
to  the  array.”  Experimental  data  is  used  with  frequencies  from  75-150  Hz  and 
ranges  from  about  3-5  km.  Thode  discusses  the  use  of  a  Radon  transform  to  fnd 
the  slopes  of  striations. 

Time-reversal  mirrors  and  adaptive  array  processing 

In  [57],  Tao  and  Krolik  used  the  waveguide  invariant  to  improve  covariance  matrix 
estimation  for  adaptive  beamforming,  in  order  to  reduce  coherent  multipath  inter¬ 
ference.  The  basic  idea  is  that  the  covariance  matrix  at  one  frequency  can  provide 
useful  information  about  the  covariance  matrix  at  another  frequency,  with  the  value 
of  (3  describing  how  the  wavenumber  differences  change  with  frequency. 

In  [37],  Kim  et  al.  use  the  waveguide  invariant  for  null-broadening  —  an 
adaptive  beamforming  technique  -  in  the  context  of  matched  f  eld  processing  in 
a  waveguide.  A  relationship  to  adaptive  time-reversal  mirrors  is  also  discussed. 
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In  [54],  Song  et  al.  use  the  waveguide  invariant  to  create  a  time-reversal  mirror 
with  variable  range  focusing.  The  idea  is  that  by  shifting  the  frequency  of  the 
received  signal  before  re-transmitting  it,  the  range  at  which  the  sound  will  focus 
will  shift  proportional  to  the  value  of  j3.  The  range  is  about  6  km  and  the  frequency 
is  about  450  Hz. 

In  [32],  Hodgkiss  et  al.  experimentally  verily  the  ideas  presented  in  [54],  and 
perform  the  TRM  focusing  at  ranges  of  up  to  30  km,  at  frequencies  around  450  Hz. 
By  comparing  the  change  in  frequency  with  the  change  in  range  of  the  focal  point, 
they  measure  /3  =  1.4. 

In  [38],  Kim  et  al.  use  the  waveguide  invariant  to  predict  the  cross-spectral 
density  matrix  (CSDM)  at  one  range  based  on  the  measured  CSDM  at  a  slightly 
different  range  from  a  ‘probe’  source  (assuming  the  depth  remains  f  xed).  This 
predicted  CSDM  is  used  to  place  a  null  with  an  adaptive  TRM. 

In  [39]  by  Kim  et  al.  changes  in  focal  depth  due  to  changes  in  frequency  and 
waveguide  depth  are  related  to  the  waveguide  invariant. 


Geoacoustic  Inversion 

In  [29],  Heaney  uses  the  waveguide  invariant  striations’  slopes  and  spacing,  along 
with  the  incoherent  transmission  loss,  for  “rapid  geoacoustic  inversion”.  In  [28], 
Heaney  extends  that  work  to  range-dependent  environments,  and  in  [30]  Heaney  et 
al.  use  a  seismic  survey  source. 


The  array  invariant  and  the  waveguide  invariant 

[43]  describes  an  invariant  property  of  acoustics  f  elds  in  a  waveguide  which  is 
dubbed  “The  array  invariant”,  which  is  used  for  passive  range  estimation  for  an 
impulsive  acoustic  source.  The  array  invariant  describes  a  property  of  the  acoustic 
signal  as  a  function  of  time  and  range  along  an  array.  The  array  invariant  differs 
from  the  waveguide  invariant  in  that  the  array  invariant  describes  the  acoustic  pres¬ 
sure  itself,  as  opposed  to  the  acoustic  intensity.  The  array  invariant  is  only  invariant 
for  modes  where  kz  has  (almost)  no  dependence  on  frequency,  as  discussed  in  [43], 
These  modes  are  usually  the  same  modes  for  which  f3  «  1. 
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A.2  Rays,  modes,  and  the  waveguide  invariant 


In  [18],  Chuprov  states  that  the  waveguide  invariant  can  be  described  using  normal 
modes  or  ray  theory.  Subsequent  research  can  be  found  in  [22]  and  [24],  the  results 
of  which  are  summarized  in  [34,  Ch.  3],  Also  see  [14]  and  [6]. 

[20,  Fig.  3.8]  provides  a  very  interesting  way  of  understanding  the  relationship 
between  rays  and  the  modal  sum,  and  that  understanding  may  be  very  useful  for  un¬ 
derstanding  how  the  observed  striations  in  I(r,  cu)  depend  on  range  and  frequency 
—  a  topic  that  is  not  well  studied  as  pointed  out  in  Sec.  2.3.3.  For  example,  (5  is 
defned  in  terms  of  a  group  of  adjacent  modes.  The  striations  in  I(r,u)  are  deter¬ 
mined  by  whichever  group  of  modes  dominates  the  modal  sum  at  particular  r  and 
oj,  as  discussed  in  Sec.  2.2.3.  Because  a  ray  represents  the  interference  between 
adjacent  modes,  ray  tracing  may  be  a  good  tool  for  determining  how  the  striations 
in  I(r,  cu)  change  with  r  and  oj.  This  author  is  not  aware  of  any  research  exploring 
this  idea  in  detail. 

The  WKB  approximation  can  be  a  “bridge”  between  normal  modes  and  ray 
theory,  as  discussed  [33,  Ch.  3]  and  even  more  so  in  [34,  Ch.  3],  [60,  Sec.  5.3]  and 
other  parts  of  [60]  also  discuss  the  relationship  between  rays  and  modes,  as  does 
[11,  Sec.  6.7.1]. 
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Appendix  B 

References  for  the  WKB 
approximation  in  ocean  acoustics 


The  WKB  approximation  is  used  in  parts  of  this  thesis  to  describe  the  horizontal 
wavenumbers  of  modes  in  ocean  acoustic  waveguides.  It  is  a  general  mathematical 
technique  for  approximating  the  solution  to  certain  types  of  differential  equations, 
and  is  discussed  in  textbooks  on  applied  math,  quantum  mechanics,  electromag¬ 
netism,  and  acoustics. 

This  appendix  is  a  far-from-exhaustive  list  of  references  that  discuss  the  WKB 
approximation  specif  cally  as  it  applies  to  ocean  acoustics. 

•  [9,  Ch.  8]  and  [8,  Ch.  23  ]  discuss  the  validity  of  the  WKB  approximation  as 
applied  to  acoustic  propagation  in  stratif  ed  media,  but  not  specif  cally  as  it 
applies  to  modes  waveguides. 

•  [8,  Sec.  48.5]  discusses  WKB  approximation  for  surface-trapped  normal 
modes  in  an  n2- linear  waveguide. 

•  [8,  Ch.  49.6]  discusses  WKB  approximation  for  surface-trapped  modes  in  an 
Epstein  waveguide. 

•  [60,  Sec.  2.9  and  5.3]  discuss  the  WKB  approximation  for  ocean  acoustic 
waveguides,  as  it  relates  to  modes  and  ray  theory. 
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•  [33,  Sec.  2.5  and  Ch.  3]  discuss  the  WKB  approximation  and  how  it  relates 
to  ray  theory.  [34,  Ch.  3]  has  even  more  material  on  how  the  WKB  approxi¬ 
mation  relates  to  ray  theory. 

•  [33,  Sec.  5.6]  uses  the  WKB  approximation  to  describe  the  mode  shapes  in 
deep-water  waveguides. 

•  [11,  Sec.  6.7]  brief  y  discusses  the  WKB  approximation,  and  also  has  a  nice 
bulleted  list  of  mode  types.  However,  it  does  not  discuss  bottom  interacting 
modes  and  thus  does  not  discuss  how  to  use  the  vertical  phase  integral  to 
obtain  the  modal  horizontal  wavenumbers.  It  does,  however,  discuss  a  good 
technique  for  normalizing  the  mode  shapes  obtained  from  the  WKB  approx¬ 
imation. 

•  [5]  discusses  how  to  use  the  WKB  approximation  to  use  calculate  the  modal 
horizontal  wavenumbers  for  modes  in  ocean  acoustic  waveguides. 

This  WKB  approximation  as  presented  in  Sec.  4.2  of  this  thesis  is  based  on  [5]. 

But  [60,  Sec.  2.9]  also  provides  most  of  the  relevant  information. 
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Appendix  C 

List  of  symbols  and  acronyms 


modal  the  adjective  for  the  noun  “normal  modes“  (e.g.,  modal  horizontal  wavenum¬ 
bers). 

JASA  The  Journal  of  the  Acoustical  Society  of  America 
MFP  matched  f  eld  processing 
SSP  sound  speed  prof  le 
TRM  time-reversal  mirror 

2D-DFT  two-dimensional  discrete  Fourier  transform 
c  speed  of  sound 
d  depth  of  ocean  waveguide 
/  scalar  acoustic  intensity 

l  mode  number  (same  as  m,  but  used  when  there  is  a  double  summation) 
m  mode  number 

n  mode  number,  but  in  terms  of  the  WKB  approximation.  See  Sec.  4.2.  In  n2- 
linear  waveguide,  n  stands  for  the  index  of  refraction. 

r  range  from  the  acoustic  point  source  in  a  waveguide 
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x  range  from  the  acoustic  line  source  in  a  waveguide.  See  Ch.  6. 

Sp  phase  slowness  (See  Sec.  2.3.3) 

Sg  group  slowness  (See  Sec.  2.3.3) 
z  depth  below  the  sea-surface 

kr  image  the  image  wavenumber  of  /win (r,  /)  in  the  r  direction.  See  Sec.  3.3.2 

kf, image  the  image  wavenumber  of  /wm(r,  /)  in  the  /  direction.  See  Sec.  3.3.2 

ft  the  waveguide  invariant  when  approximated  as  being  independent  of  the  mode 
numbers  and  frequency. 

ftmi  the  waveguide  invariant  for  a  pair  of  modes  with  mode  numbers  m  and  l. 

The  array  processing  in  Ch.  5  uses  the  following  symbols  (bold  faced  lower¬ 
case  letters  denote  vector  quantities): 

•  co  is  the  speed  of  sound  at  the  HLA  depth 

•  ko  =  is  the  wavenumber  at  the  HLA  depth 

•  J  is  the  number  of  elements  in  the  HLA,  with  each  element  numbered  as 
3  =  I ,  3  =  2,  j  =  J 

•  9l  is  the  look  direction  as  measured  from  broadside  of  the  HLA 

•  =  h)  sin  0  jj  is  the  horizontal  wavenumber  corresponding  to  the  look  di¬ 
rection 

•  v(fci)  =  [v\{ki)i  V2  (A:l),  •••  ,  uj(Ll)]t  is  the  steering  vector  correspond¬ 
ing  to  the  look  direction,  v  is  a  column  vector 

•  w  =  [w i ,  W'2 ,  •  •  •  ,  wj]T  the  array  weights  column  vector  (e.g.,  Hamming, 
Hann,  uniform) 

•  9s  is  the  bearing  of  the  acoustic  source  as  measured  from  broadside  of  the 
HLA,  in  the  plane  parallel  to  the  ocean  surface 
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•  T  is  the  wavenumber  response  function  of  the  array  processor,  which  de¬ 
pends  on  the  array  geometry  and  the  array  weights  w,  and  is  described  in 
more  detail  later. 
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16.  Abstract  (Limit:  200  words) 

The  range-frequency  waveguide  invariant  describes  striations  that  often  appear  in  plots  of  acoustic  intensity  versus 
range  and  frequency.  Signal  processing  techniques  based  on  the  waveguide  invariant  are  able  to  exploit  the  effects  of 
the  ocean  acoustic  waveguide  without  requiring  detailed  knowledge  of  the  sound  speed  profile  or  of  the  seafloor.  This 
thesis  develops  a  framework  for  understanding  and  calculating  the  waveguide  invariant,  and  uses  that  framework  to 
create  signal  processing  techniques  for  the  waveguide  invariant  which  may  be  applicable  to  passive  and  active  sonar, 
time-reversal  mirrors,  and  array  processing. 

A  method  for  passively  estimating  the  range  from  an  acoustic  source  to  a  receiver  is  developed,  and  tested  on 
experimental  data.  Heuristics  are  developed  to  estimate  the  minimum  source  bandwidth  and  minimum  horizontal 
aperture  required  for  range  estimation. 

A  semi-analytic  formula  for  the  waveguide  invariant  is  derived  using  WKB  approximation  along  with  a  normal  mode 
description  of  the  acoustic  field  in  a  range-independent  waveguide.  This  formula  illustrates  how  the  SSP  and  the 
seafloor  properties  affect  the  value  of  the  waveguide  invariant. 

Array  processing  techniques  designed  specifically  for  the  purpose  of  observing  range-frequency  striations  are 
developed  and  demonstrated.  Finally,  a  relationship  between  the  waveguide  invariant  and  wavenumber  integration  is 
derived. 
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